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Section  1.  Introduction 


The  study  of  fluid  dynamics  began  long  before  manned  flight,  and  has  progressed 
to  the  present  day  along  with  the  science  of  aerospace  engineering.  In  the  aerospace 
industry,  designs  are  pushed  to  the  limit  of  current  technological  capabilities.  The 
technology  of  flight  involves  the  understanding  of  a  complex  set  of  aerodynamic 
phenomena  that  accompany  flight.  The  complexity  of  the  aerodynamic  flows  over  flight 
vehicles  is  due  in  part  to  the  shape  and  mission  of  these  vehicles,  but  an  even  greater 
obstacle  in  fully  understanding  and  predicting  aerodynamic  behavior  is  the  complexity 
inherent  in  the  equations  governing  the  dynamics  of  fluid  flows.  Though  the 
incompressible  Navier-Stokes  equations,  which  govern  fluid  flows  at  low  Mach  numbers, 
are  simple  in  form,  the  non-linearity  in  the  equations  admits  solutions  of  a  chaotic  nature. 
It  is  the  nature  of  fluid  flows  to  break  down  and  transition  to  a  turbulent  state  at  a 
sufficiently  high  Reynolds  number.  In  a  turbulent  state,  the  full  range  of  scales  of  motion 
are  active,  from  the  length  scale  of  the  vehicle  under  consideration  down  almost  to  the 
molecular  level. 

In  1883,  Reynolds  in  his  classic  experiment  on  pipe  flows,  determined  that 
transition  to  turbulence  was  a  result  of  the  instability  of  the  smooth,  laminar  flow  with 
respect  to  small  perturbations.  Instability  can  occur  in  many  forms  and  at  many  different 
scales  in  a  complex  fluid  flow.  Boundary  layers  occur  near  vehicle  surfaces  in  regions  of 
attached  flow.  As  boundary  layers  grow  downstream  of  some  originating  location,  they 
become  unstable  and  allow  transition  to  a  turbulent  boundary  layer.  This  layer,  though 
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significantly  thicker  than  its  laminar  counterpart,  is  still  thin  with  respect  to  the 
macroscopic  features  of  the  fluid  flow.  However,  the  turbulent  boundary  layer  causes 
higher  levels  of  skin  friction  and  heat  transfer.  It  also  decreases  the  likelihood  of  flow 
separation  due  to  the  higher  level  of  kinetic  energy  in  the  boundary  layer.  Instability  is 
also  observed  in  shear  layers  separated  from  the  body  in  the  form  of  the  breakdown  of 
large  scale  structures  into  smaller,  disordered  structures.  Similarly,  vortices  are  often 
observed  to  have  a  coherent,  laminar  region  followed  by  a  breakdown  to  a  disordered 
state  that  no  longer  resembles  a  vortex.  In  all  of  these  cases,  the  breakdown  of  ordered 
flow  drastically  changes  the  nature  of  the  flow,  and  consequently,  the  aerodynamic 
performance  of  the  vehicle  under  consideration. 

In  recent  years,  great  progress  has  been  made  in  the  understanding  of  instability 
and  transition  to  turbulence.  Often  advances  are  made  by  the  investigation  of  flows 
involving  simple  geometries  rather  than  flows  over  complex  shapes.  From  the 
examination  of  transition  to  turbulence  of  these  simple  flows,  the  basic  features  of 
transition  in  more  complex  flows  can  be  studied.  In  this  study,  the  instability  and 
transition  in  channel  flows  is  investigated.  This  chapter  begins  with  the  development  of  a 
theoretical  framework  for  investigating  transition  to  turbulence  in  terms  of  the  simplified 
problem  of  instability  in  2-D  channel  flows,  presented  in  §§1.1  and  1.2.  The  role  of 
numerical  simulations  in  studying  these  flows  is  outlined  in  §1.3.  The  chapter  concludes 
with  a  summary  of  the  specific  goals  of  this  study,  in  §1.4. 

1.1.  Mechanisms  of  Transition  to  Turbulence  in  Wall-bounded 

Shear  Flows 

Though  transition  to  turbulence  in  fluid  flows  of  engineering  interest  takes  many 
different  forms  and  is  still  largely  beyond  the  realm  of  predictability,  recent  progress  has 
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been  significant.  The  early  stages  of  the  transition  process  in  certain  canonical  flows  have 
been  studied  in  great  depth  to  identify  the  physical  mechanisms  of  instability.  It  is 
expected  that  these  mechanisms  are  indicative  of  the  causes  of  transition  to  turbulence  in 
more  complex  flows. 

1.1.1.  Stages  of  Instability  and  Transition  in  Shear  Layers 

The  flows  of  interest  for  aerodynamic  applications  are  those  of  homogeneous 
incompressible  or  compressible  Newtonian  fluids.  These  flows  are  induced  by  the  motion 
of  a  body  through  the  fluid.  The  fluid,  which  must  stick  to  the  surface  of  the  body, 
undergoes  shearing  forces  in  regions  near  the  body  and  in  the  wake  of  the  body.  These 
shear  layers  are  the  sites  at  which  transition  to  a  turbulent  state  may  occur.  The  typical 
process  of  transition  in  shear  flows  can  be  considered  to  consist  of  four  stages:  namely 
receptivity,  linear  instability,  weakly  nonlinear  instability,  and  nonlinear  breakdown  (Reed 
and  Saric,  1989,  Blodgett,  1995). 

The  first  stage  of  the  transition  process,  receptivity,  is  one  of  the  least  understood. 
It  is  the  process  by  which  disturbances  become  imbedded  in  the  shear  layer.  The  source  of 
the  disturbance  may  be  the  flow  far  away  from  the  body  and  may  reach  the  body  in  terms 
of  acoustic  waves.  It  may,  for  example,  be  turbulence  in  the  freestream  flow  passing  over 
the  body,  or  it  may  be  bumps  and  imperfections  in  the  body  surface  which  provide  the 
disturbance.  The  receptivity  stage  is  complete  when  it  has  triggered  an  instability 
mechanism  in  the  shear  flow,  which  can  then  proceed  down  the  path  to  turbulence.  The 
study  of  receptivity  is  a  major  area  of  current  research.  Reviews  can  be  found  in 
Choudhari  (1993)  and  Blodgett  (1995). 

The  second  stage,  linear  instability,  is  by  far  the  most  well  understood.  Linear 
instability  is  the  study  of  the  initial  appearance  of  disturbance  waves  in  a  shear  flow. 
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These  linear  instability  waves,  though  presumed  to  be  triggered  by  some  receptivity 
mechanism,  are  greatly  determined  by  the  nature  of  the  shear  flow  itself .  One  explanation 
of  the  linear  instability  mechanism  due  to  the  equilibrium  of  forces  present  in  a  fluid.  In  a 
shear  flow,  the  fluid  moves  according  to  its  inertia,  and  its  internal  stresses  due  to  pressure 
and  friction.  These  forces  must  always  remain  in  equilibrium.  When  the  equilibrium  is 
stable,  a  small  disturbance  that  shifts  the  flow  away  from  its  equilibrium  will  not  upset  the 
balance  of  forces,  and  the  flow  will  return  to  equilibrium  after  some  time.  However,  if  the 
equilibrium  is  unstable  to  a  particular  type  of  disturbance,  it  is  possible  to  upset  the 
balance  in  such  a  way  that  the  flow  cannot  return  to  its  original  equilibrium,  even  if  the 
amplitude  of  the  disturbance  is  infinitesimally  small. 

The  particular  disturbance  to  which  the  shear  flow  is  most  unstable  is  designated  as 
the  primary  linear  instability  wave  for  that  flow.  It  is  assumed  that  the  receptivity 
mechanism  will  select  and  feed  energy  into  this  instability.  In  disturbance  environments  in 
which  many  types  of  small  disturbances  are  present,  the  most  unstable  linear  wave  will  be 
excited.  However,  it  is  possible  to  force  specific  instability  waves  in  the  flow  by  imposing 
a  disturbance  of  a  certain  type,  e.g.  of  a  certain  frequency.  In  this  case,  the  nature  of 
linear  instability  waves  are  determined  both  by  the  susceptibility  of  the  shear  flow  to 
disturbances  and  the  type  of  the  disturbances  imposed.  In  general,  though,  it  is  sufficient 
to  assume  that  the  transition  process  is  initiated  with  the  growth  of  the  most  unstable 
linear  wave.  The  search  for  the  primary  instability  of  a  shear  flow  is  performed  using 
linear  stability  theory,  which  is  outlined  in  §1.1.2. 

The  third  stage  of  the  transition  process  considers  the  evolution  of  the  primary 
instability  once  it  has  grown  to  a  noticeable  amplitude  in  the  shear  flow.  The  primary 
instability  evolves  nonlinearly  according  to  the  Navier-Stokes  equations.  It  also  interacts 
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with  any  other  disturbances  present  in  the  flow.  The  field  of  weakly-nonlinear  stability 
theory  considers  the  self-  and  mutual-interactions  of  small  amplitude  instability  waves  to 
determine  the  long-time  evolution  of  a  disturbance  once  it  is  generated.  This  approach 
will  be  discussed  in  §1.1.3. 

Once  instability  waves  grow  out  of  the  range  of  "small  amplitude,"  rapid 
breakdown  is  often  observed,  producing  the  highly  disordered  state  generically  called 
turbulence.  The  final  stage  of  transition  to  turbulence,  nonlinear  breakdown,  is  a 
combination  of  all  the  stages  following  the  weakly-nonlinear  stability  stage  until  a  fully 
turbulent  state  is  achieved.  The  only  approaches  to  understanding  this  flow  regime  are 
experimental  measurements  and  direct  simulations  of  the  Navier-Stokes  equations  using 
very  high  resolution.  Theoretical  and  statistical  tools  for  investigating  turbulence  itself  can 
be  used  to  probe  these  experimental  and  computational  measurements  of  the  transition 
regime  with  hopes  of  identifying  the  underlying  structure  of  turbulence.  From  this 
perspective,  much  work  has  been  done  in  the  investigation  of  nonlinear  breakdown.  A 
review  can  be  found  in  Kachanov  (1994). 

1.1.2.  Linear  Mechanisms 

The  earliest  method  for  examining  stability  of  fluid  flows  is  still  the  most  important 
tool  for  understanding  the  transition  process.  Linear  stability  theory  begins  by 
representing  a  solution  to  the  Navier-Stokes  equations  as  some  known  solution,  called  the 
base  state  or  base  flow,  plus  a  small  perturbation.  Neglecting  terms  that  are  nonlinear  in 
perturbation  quantities,  a  linearized  set  of  governing  equations  is  obtained.  Next,  the 
disturbance  is  assumed  to  be  periodic  in  at  least  one  direction,  typically  the  streamwise 
and/or  spanwise  directions.  The  solution  to  this  linear  system  is  then  considered  as  a 
superposition  of  Fourier  modes  in  the  periodic  direction(s).  Complex  exponential 
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variation  is  assumed  in  time  for  each  mode.  The  study  of  the  stability  of  this  linearized 
system  then  becomes  the  search  for  eigensolutions  for  each  Fourier  mode.  Each  of  the 
resulting  eigenvectors  is  a  mode  shape  of  the  perturbed  flow.  The  corresponding 
eigenvalue  contains  the  growth  rate  for  that  mode  as  its  imaginary  part,  and  the  frequency 
of  oscillation  of  that  mode  as  its  real  part.  This  modal  decomposition  of  the  flow  allows 
the  search  for  the  "inost  dangerous  mode,"  or  the  mode  with  the  most  positive  growth 
rate,  given  the  various  other  parameters  in  the  analysis,  such  as  the  Reynolds  number  and 
spatial  wavelength(s).  If  the  growth  rate  of  this  mode  is  positive,  the  base  state  is  said  to 
be  linearly  unstable.  The  standard  text  on  linear  stability  theory  is  Drazin  and  Reid  (1981). 
See  also  Bayly,  Orszag  and  Herbert  (1988)  and  Reed  and  Saric  (1989). 

Certain  pieces  of  information  are  conspicuously  missing  from  the  modal 
decomposition  explanation  of  instability.  At  first  glance,  it  seems  that  linear  theory  simply 
describes  whether  or  not  a  base  state  is  linearly  unstable;  and  if  unstable,  it  provides  the 
characteristic  shape  of  the  disturbance  that  activates  the  instability.  However,  linear 
instability  does  not  predict  the  evolution  of  the  flow,  once  departed  from  the  base  state. 
Also,  it  does  not  predict  the  response  of  the  base  state  to  a  disturbance  of  finite  size.  It 
also  is  incapable  of  estimating  the  amplitude  of  the  disturbance  that  will  be  present.  For 
these  reasons,  linear  stability  theory  alone  is  incapable  of  predicting  transition  to 
turbulence. 

However,  the  most  dangerous  mode  of  the  linear  stability  problem  is  often 
observed  in  experiments  to  dominate  all  other  modes  that  might  be  present,  and  to  grow 
to  measurably  large  amplitude  before  other  types  of  instability  can  grow  to  noticeable  size. 
This  is  evidence  that  the  most  unstable  linear  mode  plays  a  significant  role  in  the  transition 
process.  In  fact,  Henningson  and  Reddy  (1994)  show  that  linear  mechanisms  must  be 
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present  for  nonlinear  growth  of  disturbances  to  occur.  Thus,  in  spite  of  its  limitations  in 
providing  information  on  the  evolution  of  instabilities,  linear  theory  does  often  shed  light 
on  the  inception  of  the  transition  process. 

1.1.3.  Nonlinear  Growth  and  Interaction  of  Instability  Waves 

The  limitations  of  linear  theory  for  predicting  transition  led  researchers  to 
investigate  the  role  of  nonlinear  mechanisms  in  the  transition  process.  Though  many 
approaches  have  been  followed  in  considering  the  transition  process,  a  fairly  unified 
perspective  and  mathematical  treatment  of  nonlinearity  in  the  Navier-Stokes  equations  has 
become  dominant.  This  approach  is  called  weakly-nonlinear  stability  theory.  In  this 
approach,  the  primary  instability  wave  is  assumed  to  force  its  higher  harmonics  through 
the  nonlinear  convective  term  in  the  Navier-Stokes  equations.  The  theory  assumes  that 
the  primary  disturbance  is  small,  but  not  infinitesimal,  and  forces  its  harmonics  at  a  finite 
level.  The  interaction  of  the  primary  disturbance  with  its  harmonics  modulates  the  growth 
of  the  primary  as  it  feeds  some  of  its  energy  into  its  harmonics. 

Assuming  that  the  primary  disturbance  amplitude  is  very  small,  its  weakly 
nonlinear  evolution  is  given  by  the  Landau  equation  (Drazin  and  Reid,  1981).  The  Landau 
equation  is  an  ordinary  differential  evolution  equation  in  time  for  the  modulus  of  the 
amplitude  of  the  primary  disturbance,  given  its  initial  value; 

^IAP  =  2wIA|2  -  e  lAI'*  (1.1) 

where  (o  is  the  growth  rate  from  linear  theory,  and  { is  the  so-called  Landau  coefficient. 

As  the  amplitude  becomes  infinitesimal,  the  fourth-order  term  becomes  negligible,  and  the 
solution  reverts  to  that  predicted  by  linear  theory,  i.e.,  an  exponential  in  time: 
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IA|2  =  ^  lAI  =  Aoe“‘ 


(1.2) 


where  Aq  =  I  A{t  =0)  I  is  the  initial  amplitude.  For  finite  amplitude,  the  modulus  of  the 
amplitude  is  given  in  closed  form  as  a  function  of  time: 


K.  A  ^ 

- ^Aa  + 

2(0 


(1.3) 


As  shown  in  Drazin  and  Reid  (1981),  the  Landau  equation  predicts  four  categories  of 
weak  coupling  between  the  primary  mode  and  its  harmonics,  depending  on  the  sign  of  the 
(o  and  fi. 

1.  For  (o<0,  the  flow  is  linearly  stable.  That  is,  infinitesimal  disturbances  will 
decay. 

a.  For  {>0,  both  terms  on  the  right  hand  side  of  Eq.  ( 1 . 1 )  are  of  the 
same  sign,  and  a  finite  size  disturbance  will  decay. 

b.  For  f<0,  the  two  terms  of  Eq.  (1.1)  are  of  opposite  sign.  For  initial 
amplitudes  below  the  threshold  value  of: 


(1.4) 


the  disturbance  will  decay,  but  for  amplitudes  above  the  threshold 
value,  the  disturbance  will  grow  rapidly.  This  condition  is  called 
subcritical  instability.  Though  the  Reynolds  number  is  below  the 
critical  value  for  linear  instability,  the  base  flow  is  unstable  if  a  large 
enough  disturbance  is  imposed. 
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2.  For  (o>0,  the  flow  is  linearly  unstable.  Infinitesimal  disturbances  will  grow 
exponentially,  and  the  solution  will  depart  rapidly  from  the  base  flow 
solution. 

c.  For  C>0,  the  fourth-order  term  is  of  opposite  sign  to  the  linear  term. 
Therefore,  when  the  disturbance  grows  to  a  finite  amplitude,  its 
growth  will  be  slowed  until  a  threshold  value,  again  given  by 

Eq.  (1.4),  is  approached.  The  solution  will  then  remain  at  this 
solution  for  all  time.  Also,  if  the  initial  amplitude  is  larger  than  the 
threshold  value,  the  solution  will  decay  to  the  threshold  value.  This 
condition  is  called  supercritical  stability.  Though  the  flow  is 
linearly  unstable,  a  small  but  finite  amplitude  stable  solution  will  be 
obtained  when  the  base  flow  is  perturbed.  This  process  can  also  be 
considered  as  a  bifurcation  of  the  unstable  base  flow  to  a  new  flow 
in  which  the  nonlinear  interactions  are  in  equilibrium. 

d.  For  KO,  the  two  terms  of  Eq.  (1. 1)  are  of  the  same  sign.  The 
nonlinearity  will  enhance  the  growth  of  the  disturbance  as  it 
increases  in  magnitude. 

This  theory  is  the  simplest  example  of  the  theoretical  framework  used  to 
understand  nonlinear  interactions  which  occur  early  in  the  transition  process.  It  introduces 
the  concepts  of  subcritical  instability  and  supercritical  stability.  The  Landau  equation  is  a 
subset  of  a  more  general  perturbation-methods  approach  in  which  multiple  levels  of  wave 
interactions  are  considered.  Numerous  recent  studies  have  greatly  expanded  knowledge  in 
this  area.  See  Herbert  (1988b),  Reed  and  Sark  (1989),  Crouch  and  Herbert  (1993),  and 
Singer,  Erlebacher  and  Zang  (1992)  for  reviews  and  examples. 
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1.1.4.  Applications  of  Instability  Theory  in  Canonical  Flows 

How  can  the  above  concepts  be  investigated  for  the  flow  over  a  flight  vehicle  in 
real  flight  conditions?  A  simplified  understanding  of  the  stability  of  such  a  flow  can  be 
obtained  by  performing  an  experiment  on  the  flight  vehicle  to  characterize  and  measure  the 
shear  layers  on  the  vehicle.  Boundary  layers,  separated  regions,  wakes,  and  jets  are  all 
shear  layers  that  occur  over  flight  vehicles.  Is  breakdown  evident  in  these  structures?  Is 
there  a  definite  transition  location  at  which  the  flow  becomes  turbulent?  What  is  the 
spatial  extent  of  the  transition  region? 

Many  of  the  answers  to  explaining  the  observed  flow  features  lie  in  the 
investigation  of  the  stages  of  transition  in  canonical  flows  such  as  the  plane  boundary 
layer.  Receptivity  mechanisms  can  be  used  to  explain  when  and  if  instabilities  occur  in 
certain  test  conditions  and  thus  lead  to  the  hope  of  control  of  instabilities.  The  earliest 
stage  of  instability  is  often  observed  to  match  the  primary  instability  wave  predicted  by 
linear  theory.  The  subsequent  growth  of  this  wave  can  often  be  explained  by  weakly- 
nonlinear  theory.  The  nonlinear  breakdown  and  fully  turbulent  regions  can  be  understood 
in  terms  of  measurements  of  transition  and  turbulence  in  canonical  flows.  Thus  canonical 
flows  provide  answers  that  apply  to  many  types  of  more  complex  flows.  This  approach  is 
incomplete  in  explaining  all  of  the  features  of  transition,  and  is  painfully  lacking  in  the 
ability  to  predict  when  and  where  transition  will  occur  in  the  flow  over  a  complex 
geometry.  However,  it  provides  a  fundamental  understanding  of  transition,  which  is  a 
necessary  prerequisite  to  advanced  prediction  capabilities.  The  application  of  studies  of 
the  stability  of  canonical  flows  to  the  understanding  of  transition  in  complex  flows  are 
given  in  Kaiktsis,  Kamiadakis  and  Orszag  (1991)  and  Kamiadakis  and  Triantafyllou 
(1992). 


11 


Channel  flows  comprise  a  particularly  simple  class  of  canonical  flows.  Though 
channel  flows  are  internal  flows,  the  stability  and  transition  characteristics  of  channel  flows 
closely  match  those  of  boundary  layers.  Because  they  are  internal  flows,  the 
computational  and  mathematical  framework  for  their  investigation  is  much  simpler 
because  the  domain  is  small  and  wall-bounded.  These  flows  may  examined  in  greater 
depth  than  corresponding  boundary  layer  flows. 

In  this  study,  the  linear  and  weakly  non-linear  stages  in  the  growth  of  small, 

prescribed  disturbances  will  be  examined  in  the  setting  of  2-D  channel  flows.  The 

/ 

receptivity  stage  is  assumed  to  have  already  occurred,  allowing  disturbances  of  a  known 
type  to  be  present  in  the  flow.  Furthermore,  the  highly  nonlinear  stages  characterized  by 
rapid  evolution  and  breakdown  of  3-D  structures  in  the  flow  will  be  assumed  to  occur 
after  the  flow  has  evolved  according  to  linear  and  weakly-nonlinear  mechanisms.  These 
concepts  are  applied  specifically  to  channel  flows  in  the  following  section. 

1.2.  Channel  Flows 

For  the  purposes  of  this  study,  channel  flows  will  described  as  flows  between  two 
stationary  infinite  plates,  induced  by  a  mean  pressure  gradient.  The  walls  may  be  curved 
in  the  streamwise-wall  normal  plane,  but  may  not  be  curved  in  the  spanwise  direction. 

The  flow  will  be  induced  by  either  a  constant  mean  pressure  gradient  in  the  streamwise 
direction,  or  by  a  variable  mean  pressure  gradient  adjusted  to  maintain  a  constant 
streamwise  flow  rate  in  the  channel.  Flows  in  this  category  include  the  plane  channel,  or 
Poiseuille  flow,  and  curved  channel,  or  Taylor-Dean  flow,  both  of  which  are  1-D 
geometries  and  have  1-D  analytical  base  flow  solutions.  The  more  general  class  of 
periodic  2-D  channel  flows  do  not  have  analytical  base  flow  solutions.  Examples  are  plane 
and  curved  channels  with  streamwise  periodic  constrictions,  expansions,  grooves,  or  wavy 
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walls.  Examples  of  studies  of  complex  channel  flows  are  Guzman  and  Amon  (1994), 
Roberts  (1994),  Sobey  and  Drazin  (1986),  Amon  and  Patera  (1989),  Ghaddar  et  al. 

(1986),  and  Kamiadakis,  Mikic  and  Patera  (1988). 

Many  of  the  features  of  transition  observed  in  experiments  on  boundary  layers  and 
more  complex  flows  have  been  thoroughly  investigated  using  simple  channel  flow 
geometries.  Different  types  of  channel  geometries  demonstrate  very  different  routes  to 
turbulence.  Plane  channel  flow  is  one  of  the  simplest  examples  of  a  flow  which  exhibits 
subcritical  instability,  while  curved  channel  flow  exhibits  supercritical  stability.  The  initial 
nonlinear  stages  of  instability  give  plane  and  curved  channel  flows  very  different  routes  to 
turbulence.  In  fact,  these  two  canonical  flows  serve  as  examples  of  the  two  major 
categories  of  routes  to  turbulence:  the  secondary  instability  route  to  turbulence,  and  the 
Ruelle-Takens-Newhouse  scenario. 

1.2.1.  Plane  Channel 

Plane  channel  flow  is  one  of  the  canonical  settings  for  the  secondary  instability 
route  to  turbulence.  In  spite  of  the  simplicity  of  the  base  flow  solution,  the  explanation  of 
the  transition  scenario  occurring  in  plane  channel  flow  is  still  a  topic  of  debate  and 
continued  study.  However,  a  general  framework  for  the  early  stages  of  transition  in  plane 
channel  flow  has  been  agreed  upon  and  argued  to  be  typical  of  a  broad  range  of  flows 
susceptible  to  secondary  instabilities.  See  the  studies  of  Orszag  and  Patera  (1983),  Singer, 
Reed  and  Ferziger  (1989)  and  Zang  and  Hussaini  (1985,  1987,  and  1990)  and  the  review 
of  Kleiser  and  Zang  (1991). 

The  first  stage  of  the  transition  process  in  plane  channel  flow  is  the  growth  to  a 
finite  amplitude  of  a  primary  instability  mode  superimposed  on  the  1-D  parabolic  base 
flow.  Though  many  types  of  linear  instability  modes  exist  in  plane  channel  flow,  such  as 
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2-D  Tollmien-Schlichting  modes,  oblique  Tollmien-Schlichting  modes,  streamwise  center 
modes,  and  streamwise  vortex  modes,  the  base  flow  is  linearly  stable  to  all  of  these 
disturbances  at  Reynolds  numbers  at  which  transition  to  turbulence  has  been  observed  in 
experiments,  Reynolds  numbers  as  low  as  1,000.  However,  weakly  nonlinear  theory  for 
channel  flow  predicts  subcritical  instability  of  the  2-D  Tollmien-Schlichting  mode  for  an 
initial  disturbance  amplitude  above  a  threshold  value,  for  Reynolds  number  above  2,900 
and  below  the  critical  Reynolds  number  of  5,772.  Furthermore,  for  Reynolds  below 
2,900,  nonlinear  self-interaction  of  the  finite-amplitude  Tollmien-Schlichting  mode  causes 
it  to  decay  on  a  slow,  viscous  time  scale,  effectively  lowering  the  threshold  Reynolds 
number  for  which  further  bifurcations  may  occur  (Orszag  and  Patera,  1983). 

The  initial  departure  from  a  2-D  finite-amplitude  state  to  a  3-D  state  is  the  result  of 
a  linear  instability  of  the  2-D  state  to  spanwise  disturbances.  This  is  called  secondary 
instability.  Experiments  and  numerical  simulations  have  shown  that  the  spanwise, 
secondary  instability  leads  to  the  creation  of  3-D  vortical  structures  called  A  (lambda) 
vortices.  These  structures  grow  to  a  finite  amplitude,  and  then  rapidly  break  down, 
presumably  to  some  tertiary  instability.  This  process  continues  until  a  turbulent  state  is 
reached.  The  scenario  by  which  instability  modes  grow  in  time  to  a  finite  amplitude  and 
then  become  unstable  to  other  instability  modes,  with  this  process  continuing  until  full 
nonlinear  breakdown  is  achieved,  is  called  the  secondary  instability  route  to  turbulence 
(Orszag  and  Patera,  1983,  Bayly,  Orszag  and  Herbert,  1988,  Reed  and  Saric,  1989, 
Herbert,  1988b). 

The  progression  of  the  2-D  Tollmien-Schlichting  state  with  a  linearly  growing 
spanwise  instability  mode  is  not  unique,  but  is  highly  dependent  on  the  disturbance 
environment.  Nonlinear  interactions  between  various  linear  stability  modes  present  in  the 
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initial  flow  can  cause  different  3-D  states  to  result.  In  particular,  theoretical  attention  has 
been  given  to  the  interaction  of  streamwise  vortices  (linear  vortex  modes),  the  finite- 
amplitude  Tollmien-Schlichting  wave,  and  the  secondary  instability  mode  (Hall  and  Smith, 
1989,  Nayfeh  and  Al-Maaitah  1988,  Orszag  and  Patera,  1983). 

The  growth  of  the  secondary  instability  mode  is  followed  quickly,  i.e.,  on  a 
convective  time  scale,  by  breakdown  to  turbulence.  Observations  of  this  breakdown 
suggest  that  the  secondary  instability  mode  grows,  with  nonlinear  interaction  with  other 
linear  modes,  into  a  lambda  vortex  structure.  The  time  evolution  of  this  vortex  structure 
can  be  considered  in  terms  of  inviscid  vortex  filament  concepts.  The  bent  vortex  shape 
will  be  unstable  and  will  have  a  tendency  to  pull  together  like  a  hairpin,  stretching  in  the 
streamwise  direction,  and  then  breaking  apart  leaving  streamwise  streaks  at  the  wall.  This 
explanation  is  satisfying  to  those  who  observe  similar  structures  in  fully-developed 
turbulence  in  channels  and  boundary  layers. 

1.2.2.  Curved  Channel 

The  flow  in  a  channel  with  constant  radius  of  curvature  exhibits  qualitatively 
different  behavior  than  plane  channel  flow  even  for  very  slight  curvature. 

The  most  dominant  instability  mode  of  curved  channel  flow  even  with  small  levels 
of  curvature  is  a  pair  of  counter-rotating  stationary  streamwise  vortices.  This  flow,  called 
Dean  vortex  flow,  results  from  instability  of  the  base  flow  to  spanwise  disturbances.  The 
wavelength  of  the  Dean  vortex  mode  (a  single  pair  of  counter-rotating  vortices)  is 
approximately  2.5  times  the  channel  height. 

For  Reynolds  numbers  below  the  critical  value  given  by  linear  theory,  curved 
channel  flow  is  stable,  even  to  finite  size  disturbances.  Above  the  critical  Reynolds 
number,  the  base  flow  is  unstable  and  a  Dean  vortex  mode  will  grow.  However,  the 
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vortex  mode  will  saturate  at  a  small  finite  amplitude.  Thus  curved  channel  flow  follows 
the  supercritical  stability  model.  The  growth  and  final  amplitude  of  the  Dean  vortex  mode 
can  be  qualitatively  predicted  by  weakly  nonlinear  theory.  This  vortex  flow  is  stable  to 
streamwise  disturbances  for  a  large  range  of  Reynolds  number. 

The  route  to  turbulence  in  curved  channel  flow  is  currently  thought  to  follow  the 
Ruelle-Takens-Newhouse  scenario.  The  route  to  turbulence  is  understood  in  terms  of  the 
bifurcations  undergone  by  the  flow  as  the  Reynolds  number  is  slowly  increased.  At  a  fixed 
Reynolds  number  between  the  critical  Reynolds  number  and  the  value  at  which  full 
nonlinear  breakdown  is  observed,  the  flow  will  settle  into  an  intermediate  finite-amplitude 
state  dependent  on  the  Reynolds  number.  As  the  Reynolds  number  is  increased  the  flow 
will  bifurcate  from  to  a  more  complex  state.  The  Ruelle-Takens-Newhouse  scenario 
suggests  that  the  first  bifurcation  of  the  flow  is  from  the  base  state  to  a  periodic  flow  with 
frequency  Wi.  At  higher  Reynolds  number,  this  flow  will  bifurcate  to  one  with  two 
incommensurate  frequencies,  to,  and  0)2 ,  that  is  co,  is  not  an  integer  multiple  of  (£>2  and  vise 
versa.  This  results  in  a  so-called  quasi-periodic  state.  Eventually,  a  state  will  appear  with 
three  frequencies.  Dynamical  systems  theory  predicts  that  after  the  appearance  of  the 
third  incommensurate  frequency,  any  disturbance  will  trigger  chaos.  In  the  Ruelle-Takens- 
Newhouse  scenario,  a  quasi-periodic  flow  with  three  frequencies  is  rapidly  followed  by 
turbulence  (Guckenheimer,  1986). 

This  describes  the  route  to  turbulence  in  curved  channel  flow.  Starting  with  the 
Dean  vortex  flow  as  the  base  flow,  a  critical  Reynolds  number  (dependant  on  the 
curvature)  exists  above  which  the  flow  is  linearly  unstable  to  a  traveling  streamwise 
disturbance.  The  flow  at  a  particular  point  is  then  periodic  with  frequency  o)j  .  Two 
different  types  of  streamwise  instability  modes  have  been  found  (Finlay,  Keller  and 
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Ferziger,  1988,  Bland  and  Finlay,  1991).  For  lower  Reynolds  numbers,  a  long  wavelength 
disturbance  can  yield  an  "undulating"  vortex  similar  to  the  wavy  vortex  mode  found  in 
Taylor-Couette  flow.  For  higher  Reynolds  numbers  and  greater  channel  curvature,  a 
"twisting"  vortex  with  a  shorter  disturbance  wavelength  is  more  unstable  and  becomes 
dominant.  At  higher  Reynolds  number,  the  flow  in  computations  has  been  observed  to  be 
quasi-periodic  with  two,  then  three  frequencies,  suggesting  the  Ruelle-Takes-Newhouse 
scenario  for  curved  channel  flow. 

1.3.  The  Roles  of  Spatial  and  Temporal  Numerical  Simulations 

As  seen  from  the  descriptions  above,  the  transition  to  turbulence  in  wall-bounded 
shear  flows  is  a  complex  process,  one  which  has  challenged  researchers  for  over  a  century. 
Progress  in  the  qualitative  understanding  of  the  mechanisms  of  instability,  saturation, 
bifurcation  and  chaos  has  been  accompanied  by  progress  in  the  prediction  of  quantitative 
aspects  of  transitional  flow,  such  as  mode  shapes  and  growth  rates,  critical  Reynolds 
numbers  and  amplitudes  of  saturated  states.  Progress  has  been  obtained  through 
cooperation  and  dialog  between  researchers  using  many  different  techniques  for  studying 
transition.  In  many  of  the  studies  cited  above,  experimental  observations  were  used  as  the 
impetus  for  the  research.  Numerical  simulation  studies  of  transition  are  invariably 
accompanied  by  stability  theory  results  and  predictions.  Theoretical  studies  rely  heavily  on 
numerical  methods  for  solving  the  derived  equations.  All  of  the  methods  for  studying 
transition  explain  observed  behavior  in  the  language  of  mathematics  of  nonlinear  systems. 

In  this  study,  numerical  simulation  methods  are  used  to  study  the  early  stages  of 
transition  in  channel  flows.  Two  basic  approaches  may  be  taken  to  the  numerical 
simulation  of  the  stability  of  flows.  For  spatially  periodic  base  flows,  the  temporal 
simulation  method  is  used  the  most  often.  For  non-periodic  flows,  i.e.,  those  with  spatial 
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development,  the  spatial  approach  is  more  appropriate.  These  two  approaches,  along  with 
the  theoretical  framework  for  the  numerical  investigation  of  flow  instability  and  transition 
in  channel  flows  is  described  in  the  following  sections. 

1.3.1.  Assumptions  of  Theoretical  Approaches 

Theoretical  approaches  to  the  study  of  the  Navier-Stokes  equations  by  necessity 
involve  assumptions  and  approximations.  Asymptotic  methods  find  approximate  solutions 
of  the  Navier-Stokes  equations  using  a  truncated  infinite  series  in  powers  of  a  small 
parameter.  Truncation  of  the  system  results  in  simplified  equations  for  which  solutions 
can  be  sought.  Linearization  of  the  governing  equations,  and  thus  linear  stability  theory, 
results  from  the  first-order  tmncation  of  the  series.  Weakly  nonlinear  theory  is  obtained 
by  retaining  additional  terms  in  the  truncation.  These  approximations  are  valid  only  in  a 
certain  range  of  the  perturbation  parameter.  This  type  of  approximation  gives  the  bounds 
of  applicability  of  theoretical  results.  See  Crouch  and  Herbert  (1993)  for  a  complete 
example. 

Another  type  of  approximation  used  in  most  theoretical  approaches  is  the 
simplification  of  the  boundary  conditions  from  those  required  for  an  exact  solution  to  the 
problem.  In  particular,  most  stability  theories  replace  the  physical  domain  used  in 
experiments,  which  will  contain  inflow  and  outflow  regions  and  side  walls,  with  a  periodic 
domain  of  a  given  length.  That  is,  the  geometry  is  assumed  to  be  infinite  in  extent  but 
repeating  itself  in  the  streamwise  and/or  spanwise  directions  ad  infinitum.  Several 
difficulties  result  from  this  type  of  approximation.  First,  periodicity  strictly  applies  only  to 
fully  developed  flows  and  should  give  an  accurate  approximation  of  the  flow  field  only  in 
regions  far  from  the  entrance,  exit  and  side  walls.  This  is  not  a  problem  for  geometries  in 
which  the  streamwise  and  spanwise  extent  is  very  large  in  comparison  to  the  distance 
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between  channel  walls.  Second,  the  choice  of  a  periodic  wavelength  necessarily  excludes 
behavior  associated  with  other  wavelengths.  Linear  and  weakly  nonlinear  stability  theory 
considers  the  periodic  wavelength  of  the  disturbance  as  a  free  parameter.  This  allows 
selection  of  the  wavelength  in  the  solution  process,  but  it  does  not  allow  for  the  influence 
of  wavelengths  larger  than  the  selected  periodicity  to  influence  the  stability.  Third, 
periodicity  assumes  an  ordered,  repeating  structure  to  the  flow  solution.  In  most  cases 
experimental  results  show  streamwise  and  spanwise  variation  of  structures  that  are  much 
less  orderly  than  the  results  obtained  from  theory. 

A  third  type  of  approximation  used  for  theoretical  studies  is  an  assumed  form  of 
the  solution.  The  truncated  governing  equations  with  approximate  boundary  conditions 
are  solved  by  assuming  a  particular  form  of  the  solution.  The  spatial  form  of  the  solution 
is  required,  because  of  periodic  boundary  conditions,  to  be  a  wave.  The  temporal 
variation  of  the  equations  is  assumed  as  well.  Linear  stability  theory  assumes  a  complex 
exponential  variation  of  the  solution  in  time,  allowing  wave  behavior  through  the 
imaginary  part  of  the  exponent,  and  pure  growth  or  decay  through  the  real  component. 
Weakly  nonlinear  theory  results  from  modifying  this  assumption  and  creating  an  initial 
value  problem  for  the  temporal  variation.  However,  in  nonlinear  theory,  many  ad  hoc 
assumptions  must  be  made  as  to  the  relative  contribution  of  various  terms  in  determining 
the  temporal  evolution  of  the  perturbation. 

The  above  approximations  are  the  art  and  essence  of  using  theoretical  approaches; 
the  computed  solution  is  only  as  good  as  the  approximations  used  in  obtaining  it. 
However,  these  approximations  allow  the  standard  mathematical  tools  of  Fourier  analysis, 
eigenvalue  methods,  and  nonlinear  dynamical  systems  analysis  to  be  applied  to  the 
incompressible  Navier-Stokes  equations.  They  also  allow  parametric  variation  and  study 
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of  general  classes  of  flow  solutions.  Most  of  the  advances  in  the  understanding  of 
instability  and  transition  have  resulted  from  the  use  of  such  methods  in  conjunction  with 
experiments  or  numerical  simulations. 

1.3.2.  Numerical  Simulation  of  Spatially  Developing  Flows 

The  most  natural  approach  for  simulating  an  experimentally  observed  flow  using  a 
computational  method  is  to  attempt  to  match  the  experimental  setup.  Inflow  and  outflow 
regions  are  required,  as  well  as  side  walls.  The  inflow  condition  must  match  as  closely  as 
possible  the  inflow  in  the  experiment;  the  incoming  flow  will  have  certain  velocity  profile 
and  disturbance  characteristics.  The  outflow  conditions  must  also  be  selected  so  that 
numerical  difficulties  are  not  encountered  which  would  cause  different  behavior  than  that 
observed  in  the  experiment.  Typically,  an  idealized  version  of  the  experimental  boundary 
conditions  is  used.  The  inlet  profile  is  usually  assumed  to  be  a  fully-developed  channel 
flow,  plus  a  disturbance.  The  disturbance  is  usually  a  traveling  wave  solution  obtained 
from  the  linearized  Navier-Stokes  equations,  or  a  combination  of  such  solutions.  Random 
disturbances  have  also  been  used.  The  outflow  condition  usually  requires  special 
treatment,  such  as  modification  of  the  governing  equations  in  a  buffer  domain  near  the 
outflow,  so  that  disturbances  will  travel  smoothly  out  of  the  domain  without  reflection 
(Blodgett,  1995,  Danabasoglu  and  Biringen,  1991). 

The  methodology  for  simulating  instability  in  such  flows  is  to  first  obtain  a  base 
flow,  and  then  to  perturb  that  base  flow  to  determine  its  stability.  The  base  flow  is 
obtained  by  performing  the  simulation  without  the  disturbance.  A  time-asymptotic  state, 
usually  steady,  is  obtained  from  the  numerical  simulation  starting  from  some  suitable  initial 
condition.  Then  a  second  simulation  is  performed  with  the  added  disturbance.  The  effect 
of  the  disturbance  on  the  base  flow  is  obtained  by  subtracting  the  first  solution  from  the 
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second.  The  problem  has  the  character  of  a  boundary  value  problem  for  the  disturbance, 
in  that  the  disturbance  is  generated  through  an  inflow  boundary  condition  and  is  allowed 
to  grow  or  decay  spatially  in  the  domain. 

Using  the  spatial  simulation  method,  the  disturbances  will  be  observed  to  grow  or 
decay  as  they  travel  from  the  inflow  to  the  outflow.  If  the  base  flow  is  unstable  to  the 
perturbation,  the  onset  of  three-dimensionality,  quasi-periodicity,  etc.,  will  be  observed  in 
the  spatial  development  of  the  flow.  This  is  in  marked  contrast  to  the  theoretical 
approach,  in  which  the  entire  flow  bifurcates  to  a  new  state  when  the  previous  state 
becomes  unstable  due  to  the  global  nature  of  the  instability  growth  assumed  through  the 
periodicity  condition.  However,  the  spatial  development  can  be  related  to  the  temporal 
bifurcation  by  realizing  that  the  temporal  approach  describes  the  evolution  of  a  single 
wave  traveling  downstream  at  the  convective  wave  speed  of  the  flow. 

Though  the  spatial  approach  seems  ideal,  it  is  often  impractical  due  to  the 
requirement  that  the  domain  must  be  long  enough  to  contain  all  of  the  desired 
development  of  the  disturbance  from  its  inflow  value  to  its  finite  growth  and  possible 
bifurcation.  However,  certain  features  of  the  flow  can  be  observed  in  spatial  simulations 
that  are  not  adequately  described  by  the  temporal  approach.  Solutions  can  exhibit  less- 
ordered  aperiodic  structure  observed  in  experiments.  An  example  of  this  is  the  study  of 
Guo  and  Finlay  (1994)  on  spatially  developing  curved  channel  flow,  in  which  Dean 
vortices  are  observed  to  slowly  merge  and  split  as  they  travel  downstream,  as  seen  in  the 
experiments  of  Ligrani  (1992),  Finlay  and  Nandakumar  (1990)  and  Ligrani  and  Niver 
(1988).  Thus  the  Dean  vortices  are  not  strictly  parallel  to  one  another  as  is  assumed  in  the 
temporal  approach.  These  issues  are  addressed  in  Kleiser  and  Zang  (1991),  Blodgett 
( 1 995)  and  Danabasoglu  and  Biringen  (1991). 
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1.3.3.  Temporally  Developing  Flows 

Though  spatial  simulations  of  channel  flow  instability  are  necessary  for  studying 
certain  aspects  of  the  transition  process,  temporal  simulations  are  predominantly  used.  In 
the  temporal  numerical  simulation  method,  the  study  proceeds  much  more  closely  with 
theory.  Since  theoretical  approaches  often  assume  periodic  boundaries  in  space,  the 
temporal  numerical  simulations  can  be  used  to  provide  a  fully  nonlinear  solution  for 
comparison  with  weakly  nonlinear  theories.  The  validity  of  the  approximations  (other 
than  simplified  boundary  conditions)  made  in  the  theoretical  approach  are  then  revealed. 
Thus  temporal  numerical  simulation  is  used  more  often  than  spatial  simulation  when 
investigating  the  nature  of  bifurcations  and  the  transition  scenario. 

The  temporal  numerical  simulation  approach  is  quite  different  than  that  of  the 
spatial  approach  due  to  the  different  role  of  the  boundary  and  initial  conditions.  In  the 
temporal  approach,  the  disturbance  is  included  with  the  base  flow  as  the  initial  condition. 
The  periodic  boundary  conditions  insure  that  when  this  disturbance  travels  out  of  the 
domain,  it  re-enters  through  the  inflow.  The  base  flow  is  obtained  as  a  time-asymptotic 
state  with  the  same  boundary  conditions,  with  no  disturbance.  In  the  case  of  plane  or 
curved  channel  flow,  a  fully-developed  profile,  either  analytical  or  numerically 
approximated,  may  be  used  as  the  base  flow.  The  perturbation,  as  well  as  the  periodic 
wavelength,  is  obtained  from  the  linearized  Navier-Stokes  equations.  The  temporal 
evolution  of  the  disturbance  is  then  investigated.  If  the  base  flow  is  unstable  to  the 
disturbance,  the  entire  flow  will  bifurcate  to  a  new  stable  flow  or  will  breakdown  in  time 
and  asymptote  to  a  fully  turbulent  state.  The  solution  method  has  the  character  of  an 
initial  value  problem,  with  all  of  the  dynamics  of  the  flow  determined  by  the  nature  of  the 
initial  conditions.  See  Kdeiser  and  Zang  (1991)  for  a  thorough  discussion. 
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The  temporal  approach  has  several  significant  advantages  over  using  the  spatial 
approach.  The  solution  may  be  obtained  with  much  less  computational  resources  since  the 
domain  length  is  one  wavelength  of  the  disturbance.  Therefore,  more  advanced  stages  of 
the  transition  process  can  be  studied  for  flows,  such  as  plane  channel  flow,  in  which  fully 
turbulent  flow  will  eventually  occur.  Another  significant  advantage  is  that  the  stability  and 
bifurcation  behavior  can  be  thoroughly  studied  by  varying  the  Reynolds  number  of  the 
flow.  The  critical  Reynolds  number  of  the  base  flow  can  be  determined  by  slowly 
increasing  the  Reynolds  number  until  growth  of  the  disturbance  is  observed.  Subsequent 
bifurcations  may  be  obtained  by  further  increasing  the  Reynolds  number.  In  this  manner 
the  critical  Reynolds  numbers  and  nature  of  the  bifurcations  can  be  observed.  See  Singer, 
Erlebacher  and  Zang  (1992)  and  Guzman  and  Amon  (1994)  for  examples. 

In  summary,  the  mathematical  tools  used  for  understanding  transition  are  usually 
framed  in  terms  of  temporal  growth  of  disturbances  in  spatially  periodic  flows.  Temporal 
simulation,  therefore,  closely  mirrors  the  theoretical  approach  and  places  the  results  in  the 
framework  of  bifurcation  theory,  or  the  theory  of  nonlinear  dynamical  systems. 

1.4.  Goals  of  this  Study 

As  was  discussed  above,  numerical  simulation  of  the  Navier-Stokes  equations  with 
periodic  boundary  conditions  can  be  effectively  used  for  investigating  the  stages  of 
transition  in  general  2-D  channel  flows.  It  was  also  observed  that  the  transition  scenario 
for  channel  flows  has  application  to  understanding  transition  to  turbulence  in  shear  flows 
in  general  and  to  engineering  applications  in  aerodynamics.  The  long  term  purpose  of  this 
research  is  to  investigate  the  routes  to  turbulence  in  channel  flows  with  the  eventual 
application  of  explaining  transition  phenomena  in  more  general  flow  configurations. 
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The  specific  purpose  of  this  study  is  to  develop  a  tool  for  investigating  the  stability 
of  the  incompressible  flow  in  2-D  channels  using  a  3-D  numerical  simulation.  In  Section 
2,  the  formulation  of  the  equations  governing  flows  in  2-D  channels  will  be  described, 
including  the  3-D  Navier-Stokes  equations  in  general  curvilinear  coordinates,  and  the  2-D 
linearized  Navier-Stokes  equations  for  plane  and  curved  channel  flows,  for  use  in  linear 
stability  theory.  In  Section  3,  a  spectral  numerical  simulation  method  will  be  described  for 
solving  the  3-D  Navier-Stokes  equations.  The  methods  for  using  the  results  of  linear 
stability  theory  to  obtain  the  initial  condition,  and  the  extraction  of  stability  information 
from  the  computed  results  are  also  discussed.  In  Section  4,  the  development  and 
validation  of  a  computer  code  for  simulating  2-D  flow  in  plane  and  curved  channels  is 
described  along  with  the  flow  results  for  these  two  geometries.  The  results  primarily 
validate  the  simulation  methodology,  but  they  also  show  some  of  the  interesting  effects  of 
nonlinearity  in  the  early  stages  of  transition.  Particularly,  a  finite-amplitude  equilibrium 
state  is  computed  for  plane  channel  flow  at  a  subcritical  Reynolds  number.  Section  5 
contains  conclusions  and  recommendations  for  further  study. 


Section  2.  Mathematical  Formulation 


The  incompressible  Navier-Stokes  equations  describe  the  flow  of  a  liquid  or  a  low- 
speed  gas  which  is  homogeneous  and  Newtonian.  Also,  it  is  assumed  that  viscosity  is 
constant,  or  equivalently,  that  temperature  gradients  in  the  flow  are  negligible  and  do  not 
need  to  be  computed  or  considered  to  affect  the  viscosity.  The  incompressible  Navier- 
Stokes  equations  consist  of  the  conservation  of  linear  momentum,  which  is  the  application 
of  Newton's  second  law  to  an  infinitesimal  fluid  volume,  and  the  conservation  of  mass,  or 
the  continuity  equation. 

The  Navier-Stokes  equations  are  formulated  in  terms  of  general  curvilinear 
coordinates,  with  application  to  a  general  2-D  channel,  in  §2. 1 .  The  simplified  forms  of 
the  Navier-Stokes  equations  for  plane  and  curved  channel  flow  are  also  given.  In  §2.2, 
the  linear  stability  equations  are  derived  for  the  plane  and  curved  channel  by  linearizing  the 
Navier-Stokes  equations  about  a  base  state  and  applying  a  Fourier  mode  disturbance.  The 
numerical  aspects  of  solving  these  equations  and  using  them  for  the  study  of  instabilities  in 
channel  flows  is  reserved  for  Section  3. 

2.1.  Governing  Equations  and  Base  Flow  Solutions 

The  Navier-Stokes  equations  in  vector  form  are  deceptively  simple.  This  section 
applies  these  equations  to  a  general  2-D  channel  geometry,  which  requires  a  general 
curvilinear  coordinate  transformation  in  the  plane  of  the  streamwise  and  wall-normal 
directions.  The  channel  is  assumed  to  be  infinite  with  no  curvature  in  the  spanwise 
direction,  thus  allowing  a  simple  linear  transformation.  The  purpose  of  the 
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transformations  is  to  map  the  geometry  onto  a  computational  box  which  can  be  used  to 
solve  the  flow  in  any  2-D  channel  geometry.  The  computational  domain  will  be  further 
discussed  in  the  numerical  formulation  of  Section  3. 
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In  addition,  it  is  necessary  to  derive  the  Navier-Stokes  equations  for  the  specific 
cases  of  the  plane  channel,  and  the  constant  curvature  channel  (called  a  curved  channel), 
so  that  the  linearized  stability  equations  for  these  cases  can  be  derived  in  §2.2.  These  two 
special  cases  have  analytical  solutions  which  are  used  in  the  stability  analysis  and  in  the 
initial  condition  for  the  numerical  simulation.  These  solutions  are  also  derived  in  this 
section. 

2.1.1.  Vector  Form  of  Incompressible  Time-Dependent  Navier-Stokes  Equations 
The  non-dimensional  time-dependent  Navier-Stokes  equations  governing  the 
incompressible  flow  of  a  homogeneous  fluid  are: 


—u+^-iuu)  =  -Vp+— V^«+0 
dt  Re 

V-  M  =  0 


(2.1) 


where,  lengths  are  non-dimensionalized  by  an  appropriate  reference  value  ,  the  velocity, 
M,  is  non-dimensionalized  by  the  reference  velocity  ,  and  the  pressure  ,  p  ,  is  non- 
dimensionalized  by  U^.  The  density  of  the  fluid  is  constant  and  has  been  non- 
dimensionalized  by  the  reference  value  so  that  p  =  1 .  The  Reynolds  number.  Re,  is 
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where  |a  is  the  dimensional  viscosity  of  the  fluid.  The  vector  Q  is  reserved  for  use  as  a 
source  term  for  forcing  the  solution  during  the  validation  phase  of  the  study. 

The  domain  of  interest  is  the  region  between  two  walls  infinite  in  span.  The 
domain  is  defined  in  terms  of  the  computational  coordinates  1=1, 2, 3. 

•  Streamwise  direction;  e  [0, 27t] 

•  Spanwise  direction:  E  [0,  In] 

•  Wall-normal  direction:  e  [-1,  1] 

The  boundary  conditions  for  channel  flow  problems  are  periodic  in  two  directions: 

u(V  -1)  =  u 

.S  .  a;  “towr  wall 

=  U  ,, 

^  upper  wall 

Note  that  pressure  is  required  to  have  a  net  pressure  gradient  in  the  streamwise  direction 
to  allow  a  net  flux  of  mass  in  the  streamwise  direction  to  occur.  Therefore,  a  known 
constant  pressure  gradient  is  imposed,  and  the  pressure  is  split  into  a  known  non-periodic 
component,  P ,  which  forces  the  flow  in  the  streamwise  direction,  and  a  component,  p  , 
which  is  periodic  in  ^  ^ .  Note  that  P  is  chosen  so  that  its  gradient  is  constant  in  ^  .  Also 
note  that  no  wall  boundary  conditions  are  needed  for  pressure.  Finally,  since  the  pressure 
appears  only  in  terms  of  its  gradient,  an  additional,  arbitrary,  condition  is  required.  In  this 
study,  the  condition  that  is  specified  on  pressure  is  that  the  mean  pressure  in  the  domain  is 
zero.  These  conditions  are  expressed  mathematically  as: 
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p  ^  P  +  p 

V  P  •  e ,  =  const 

‘  (2.3) 

II L 

2.1.2.  Generalized  Curvilinear  Coordinates 

A  generalized  two-dimensional  transformation  is  used  to  map  the  physical  domain, 
given  by  Cartesian  coordinates  x ' ,  into  the  computational  domain,  given  by  computational 
coordinates  ^ ' .  Use  is  made  of  an  intermediate  coordinate  system  given  by  r|',  which  is 
related  linearly  to  the  computational  coordinate  system.  The  domain  of  each  coordinate 
is: 

^  f^r  /={  1,2,3} 

’  llinin+V  *={1’2,3} 

(2.4) 

,  27t] 

,  271] 

,  1] 

The  first  transformation,  ,  relates  a:  ‘  to  r| '  and  is  given  analytically  by: 

r|'  =  ri*(A:\A:^) 

Tr-  '  n'  -  TiLn  =  ^  ^  (2.5) 
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where  the  spanwise  computational  coordinate  is  an  arbitrary  linear  scaling  of  the  Cartesian 
coordinate  .  The  linear  transformation  Tj  then  relates  r| '  to  ^  ‘ : 


irl  271  ,  1  1- 

nl  271:  .  2  2  X 

^  =  —  (  11  -  Tlmin) 

f  ( *1’  -  nL) 


(2.6) 


Covariant  and  contravariant  base  vectors  are  defined,  respectively,  in  terms  of  the  above 
transformations: 


= 


dx^ 


e  = 


dx-' 


^J 


(2.7) 


where  represents  the  Cartesian  unit  basis  { ij,k }  for  i={  1,2,3},  respectively.  The 
Cartesian  base  vectors  can  also  be  expressed  in  terms  of  covariant  and  contravariant  base 
vectors: 


(2.8) 


Summation  notation  will  be  used  consistently  in  this  document  unless  otherwise  noted. 
The  covariant  and  contravariant  metric  tensors  are  defined  in  terms  of  the  base  vectors: 
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gij  =  e.;ej  g‘j  =  e‘-ej 

=  =  g.J  =  (2.9) 

matrix  {g..}  =  [matrix  {g'^}]'‘ 


The  Jacobian  of  the  transformation, 

j  _  d{x^,x^yK^) 


(2.10) 


is  related  to  the  determinate  of  the  co variant  metric  tensor  g  =  detlgy }  =  .  The  scaling 

of  the  volume  of  curvlinear  space  relative  to  Cartesian  space  is  given  by  \/g  =  1 7 1 . 

A  vector  field  such  as  the  velocity  field  has  several  representations  in 

the  general  coordinate  system: 

u  =  u‘  e.  Cartesian 

u  -  u'  e.  Contravariant  (2.11) 

u  =  u.  e'  Covariant 

Christoffel  symbols  arise  when  derivatives  of  base  vectors  are  required: 
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d  dx^  ^  dx^ 


where  =  g‘'" 


d^‘d^j  d^‘ 
d^x'^  dx'^ 


r^e 

ij  m 


ae'a^'  a^' 


a^  . 

'^k 


a^‘ 


n 


-  a;?, 

a^' 


where  AJ  = 

a^ 


(2.12) 


where  is  a  Christoffel  symbol  of  the  second  kind,  and  is  symmetric  with  respect  to  i 
andy.  Second  derivatives  of  base  vectors  are  also  required.  The  symbols  are  also 
symmetric  with  respect  to  i  and  j. 

2.1.3.  Contravariant  Form  of  Navier-Stokes  Equations 

To  apply  the  generalized  coordinate  system  to  the  governing  equations  Eq.  (2.1), 
the  following  differential  operators  are  needed: 


^  V  ^ 

V- 

V  VX; 


AdS 


'  r ' 

c'- 

e‘x 


—A 

ae' 


conservative  form 


non -conservative  form 


(2.13) 


where  in  the  case  of  the  gradient,  A  is  an  arbitrary  scalar,  vector,  or  tensor,  and  in  the 
divergence  and  curl  is  an  arbitrary  vector  or  tensor.  Applying  the  general  coordinate 
transformation  to  Eq.  (2.1)  along  with  the  pressure  splitting  from  Eq.  (2.3)  yields: 
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-m  d 


-e 


-p-e 


dP  1  1  a 


a? 


Re 


v^a^' 


'/gg 


ijA 


-U 


+  Q  (2.14) 


for  the  momentum  equation.  In  order  to  place  the  vector  momentum  equation  into  a  form 
in  which  it  can  be  solved,  it  must  be  separated  into  a  system  of  scalar  equations.  The 
vector  u  must  be  replaced  by  its  contravariant  components,  and  the  derivatives  of  base 
vectors  must  be  expanded  in  terms  of  Christoffel  symbols.  The  contravariant  components 
of  the  momentum  equation  become: 


Momentum  -  e ' 


-u 


1 

( 

1  a 

Re 

fgdf  '  "  "  a^'  "  ae' 


(2.15) 


^jk 


'■ijk 


+  12' 


for  m=  1,2,3.  The  continuity  equation  in  terms  of  contravariant  velocities  is: 


J_A 

fgd^‘ 


( v/gM '  )  =  0 


(2.16) 


2.1.4.  Plane  Channel  Geometry 

The  plane  channel  geometry  is  obtained  using  a  Cartesian  coordinate  system,  that 
is,  the  intermediate  coordinate  system  is  identified  with  Cartesian  coordinates,  riW  for 
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/={  1,2,3}.  The  Cartesian  coordinates  are  non-dimensionalized  by  the  channel  half-height 
h.  Thus  the  transformation  between  computational  coordinates  and  Cartesian  coordinates 
is  given  by: 


T.2  :  ^ 


1 1 

;c 

.  2^  ^2 


1^-2  Jc^  +  1 


(2.17) 


for  a  Cartesian  domain: 

,  L  .] 

x'^g[Q  ,  L  2]  (2.18) 

x^G[-l  ,  1] 


The  reference  velocity  for  the  plane  channel  flow  is  taken  to  be  the  maximum  base  flow 
velocity,  Thus  the  Reynolds  number  is  Re  =  hl\y . 

The  metrics  are  easily  determined  from  the  above  transformation.  All  of  the 
metrics  and  base  vectors  are  constant  and  the  derivatives  of  base  vectors  are  zero. 
Substituting  these  metrics  into  Eqs.  (2.14)  and  (2.16)  yields  the  contravariant  momentum 
and  continuity  equations  for  the  plane  channel  geometry  in  terms  of  the  generalized 
coordinates  ^  ‘ .  Also  of  interest  are  the  governing  equation  formulated  in  terms  of  the 
original  Cartesian  coordinates  x' .  These  equations  are  used  to  determine  analytical 
solutions  of  the  equations  and  to  formulate  the  linear  stability  equations  in  the  next 
section.  The  Cartesian  governing  equations  are: 
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^  ^  tn  .  ^  ^  i  ^  m 

- m'”+ - m'm"’  =  -- 


dt  dx  ‘ 


dx  ‘ 


-p-8 


dP  1  d  d 


Im 


dx  *  P^  dx‘  dx 


-u'"  +  Q‘ 


5  -  /•  n 

—  u  =  0 


(2.19) 


dx‘ 


The  above  system  of  equations  with  domain  given  by  Eq.  (2.18)  has  an  analytical  solution 
for  flow  in  the  streamwise  direction  with  no-slip  walls: 


u  =  U(x^)e,  =  (l-(x^)^)e^ 


dP 

dx^ 


Re 


(2.20) 


2.1.5.  Curved  Channel  Geometry 

To  obtain  the  curved  channel  geometry,  a  polar  coordinate  system  is  used.  The 
intermediate  coordinate  system  t) '  is  identified  with  polar  coordinates, 
{ri‘=0,r|^=y,ri^=r} ,  with  domain 


0e[e  .  ,  0  .  +Lft] 

^  min  ’  min  0  J 


y^[o  ,  T  ] 

rG[i?^-l  ,  R^  +  l] 


(2.21) 


where  is  the  mean  radius  of  curvature  of  the  channel.  Again  all  lengths  are  non- 
dimensionalized  by  the  channel  half-height  h.  The  standard  transformation  between  polar 
coordinates  and  Cartesian  coordinates  is: 
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' 

3) 

0  =  tan"^ 

X 

x^  =  r  cos(0) 

9 

►  r,'*  =  ^ 

=  y 

y  =  X- 

x^  =  r  sin(0) 

+(x^)^ 

(2.22) 


along  with  its  inverse.  The  transformation  between  computational  coordinates  and  polar 
coordinates  is  a  linear  transformation: 


=  ^(0-0  .  ) 

s2  27t 

y 

=  r-R, 


(2.23) 


The  metrics  relating  computational  coordinates  to  Cartesian  coordinates  are  determined 
from  transformations  Tj  and  T2.  Note  that  by  using  the  standard  polar  coordinate  system 
the  direction  of  increasing  ^ '  is  from  right  to  left  (for  0  in  the  first  two  quadrants).  For  the 
results  shown  in  this  study,  the  orientation  of  the  plots  is  reversed  so  that  ^ '  increases 
from  left  to  right. 

In  addition  to  the  governing  equations  in  terms  of  computational  coordinates  and 
contravariant  velocities,  the  equations  for  curved  channel  flow  in  polar  coordinates  are 
needed.  For  this  purpose,  an  alternate  representation  of  the  velocity  vector  is  used: 


35 


u  =  u^Cq  +  +  u’'e^ 

where  =  ^2/1^21 

~  ^3  ^  I  ^3  I 

therefore 

M  '■  = 


(2.24) 


These  are  the  polar  velocity  components  typically  used  in  the  literature  and  are  simply  a 
scaling  of  the  contravariant  components.  The  polar  momentum  equations  are  obtained  by 
taking  the  dot  product  of  the  vector  momentum  equation,  Eq.  (2.1),  with  each  of  the  polar 
unit  base  vectors  {cg,  .  The  resulting  momentum  equations  are 


d  0 

— u  +u 
dt  dr 


3m ®  du^ 


1 


r  30 

,0 


+  M- 


,3m 


®_m®m'‘  _  _l  dp  _  I  dP  ^ 
dy  r  r  dQ  r  dd 

32., 0  ;32„0  o  :3.,  /•  ..el 


3"m«  13m«  1  3^m”  3^m”  2  3m'' 

- + - + - + - + - 


M 


Re 

[  dr^ 

CD 

r^  302 

3y2 

r2 

30 

r2 

c 

>  r  , 

-1J  +1J 

.3m  '' 

+ 

du'' 

+  7J 

ydu'' 

Me2 

. 

3 

I4r  '  14- 

t 

dr 

1  f/l 

r  30 

dy 

r 

dr 

1 

d^u 

.  •+• 

1  3m'' 

1  d^u  ''  d^u  '' 

-1-  - 

u 

_  2 

3m® 

Re 

dr^ 

r  dr 

r2  302 

dy^ 

'r2 

30 

+  e' 


+  <2'' 


3  y  du^  ydu^  dp 

— M^  +  m'^ - + - +M^ -  =  -— + 

dt  3r  r  30  3y  3y 


+  Q 


1 

d^u^ 

Re 

dr^ 

r  dr  Qy 


(2.25) 


and  the  continuity  equation  is 
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\  du^  \  d  I  r  \  du^ 

— —  + - (  ru  )  + -  =  0 

r  dd  r  dr  dy 


(2.26) 


The  above  system  of  equations  has  an  analytical  solution  for  flow  in  the  streamwise 
direction  with  no-slip  walls: 


u  =  U^{r)eQ  =  Dr  Inr  +  C  +  E-^ 


where 


^ 

dQ  Re 

dP  ^  {U^Y 

dr  r 


C  =  — —  Inti  +  In/? 
4R.  ' 


2 

E  =  -Rf-^  Inri 
'  4R 

C 


/?,  =  R^-1 


K  =  K  +  ^ 

o  c 


(2.27) 


R^  +  £'lnri 

See  Singer,  Erlebacher  and  Zang  (1992).  To  be  consistent  with  the  literature,  the 
reference  velocity  for  curved  channel  flow  is  taken  to  be  the  mean  base  flow  velocity. 


f  U\r)dr 


(2.28) 


Thus  the  Reynolds  number  is  /?e  =  . 
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2.1.6.  Forcing  Function  Analysis 

The  purpose  of  the  forcing  function  Q  is  described  in  this  section.  Since  few 
analytical  solutions  of  the  Navier  Stokes  equations  are  known,  it  is  helpful  for  purposes  of 
code  validation  to  provide  an  analytical  solution  and  require  that  it  be  a  solution  to  the 
governing  equations  by  choosing  a  suitable  forcing  function.  In  general,  this  can  be  done 
by  prescribing  the  velocity  vector  u=U,^^,  in  any  coordinate  system,  along  with  the  pressure 
p  =  In  this  study,  only  divergence-free  test  solutions  I/, were  selected.  After 
substituting  the  test  solution  into  the  momentum  equation,  the  forcing  function  vector  Q 
is  chosen  to  be  the  residual  of  the  equation.  In  vector  form,  this  can  be  written  as: 


Q  ^  —u  +  V-(m«)  +Vp  -  —V^u 
dt  Re 


(2.29) 


In  this  study,  the  analytical  test  solution  was  substituted  into  either  the  Cartesian 
momentum  equations,  Eq.  (2.19),  or  the  polar  momentum  equations,  Eq.  (2.25).  The 
forcing  function  Q  was  evaluated  symbolically  using  Mathematica  software  (Wolfram, 
1991)  and  stored  within  the  code.  During  the  simulation,  the  forcing  function  vector  Q  ,  is 
a  known  function  of  space  and  time.  Adding  this  term  does  not  change  the  nature  of  the 
Navier-Stokes  equations  as  a  system  of  differential  equations. 

2.2.  Linear  Stability  Theory 

The  linear  stability  equations  for  plane  and  curved  channel  flows  are  obtained  using 
both  the  Navier-Stokes  equations  linearized  about  some  known  base  flow,  and  a  given 
form  for  disturbances  about  that  base  flow.  An  abbreviated  description  of  linear  stability 
theory  for  plane  and  channel  flows,  along  with  the  derivation  of  the  stability  equations  and 
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the  Orr-Sommerfeld  equation  are  described  below.  For  more  complete  details,  refer  to 
Drazin  and  Reid  (1981). 

2.2.1.  General  Parallel  Flow  Formulation 

Linear  stability  analysis  is  based  upon  a  series  of  assumptions  and  approximations 
about  perturbations  to  some  known  solution  of  the  governing  equations.  The  known 
solution  is  the  base  flow  or  basic  state,  u=Ui,  and  p=Pi, .  The  stability  of  this  base  flow  to 
periodic  wave  disturbances  in  a  given  direction  ^ or  ^  ,  is  investigated.  In  this  study, 

^  may  be  either  ^ ,  or  ^  ^ ,  but  not  ^  ^  since  the  walls  at  ^  ^  - 1  and  1  prohibit  a  periodic 

solution.  The  base  flow  is  assumed  to  be  parallel  in  the  wave  direction,  that  is 
d/d^Uf^  =  dld^Pf^  =  0  in  the  derivation  of  the  following  stability  equations.  Given  a 
basic  state  and  wave  direction,  the  following  steps  lead  to  linear  stability  results. 

First,  a  perturbation  is  applied  to  the  basic  state.  For  an  infinitely  small 
perturbation,  the  temporal  variation  may  be  assumed  to  be  complex  exponential  in  time 
with  both  the  velocity  vector  and  pressure  having  the  same  temporal  variation.  The 
perturbation  is  assumed  to  be  periodic  in  the  ^"'-direction  with  wavelength  X  =  2%! a 
where  a  is  the  wavenumber,  assumed  to  be  real  in  this  study.  This  puts  the  perturbation  in 
the  form: 


E 

p'-  t 


u  =  Uf^+u' 

P  =  Pb*p' 

Real  I  }  (2.30) 

Real  { } 
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where  and  are  the  complex  mode  shapes.  The  above  substitutions  are 

made  into  the  governing  equations.  Since  the  base  flow  is  known  to  be  a  solution  to  the 
governing  equations,  the  governing  equations  with  the  base  flow  substituted  are 
subtracted  from  the  perturbed  equation.  The  equations  are  then  linearized,  i.e.,  any  term 
quadratic  in  perturbation  quantities  (which  are  assumed  infinitely  small)  is  neglected. 

The  linear  stability  equations  are  obtained  analytically  by  taking  derivatives  with 
respect  to  time  and  Since  the  equations  are  linear,  each  mode  of  the  Fourier  series  may 
be  considered  independently.  Typically,  the  Fourier  mode  number,  k,  is  dropped  in  the 
development  of  the  linear  stability  equations.  Actually,  it  is  just  absorbed  into  a.  The 
system  of  linear  differential  equations  with  prescribed  boundary  conditions  only  has  a 
solution  for  certain  combinations  of  the  various  parameters  in  the  equations,  namely  the 
Reynolds  number  Re,  spatial  wavenumber  a  (actually  ka)  and  complex  temporal  growth 
rate  to.  In  this  study,  the  Reynolds  number  and  spatial  wavenumber  are  prescribed.  The 
unknown  parameter,  the  complex  temporal  growth  rate,  can  then  be  determined  as  an 
eigenvalue  of  the  system  of  differential  equations.  The  eigensolution  corresponding  to  a 
particular  eigenvalue  consists  of  the  mode  shapes  u  and  p .  The  linear  stability  equations 
for  streamwise  disturbances  in  plane  and  curved  channel  flow  are  given  in  the  following 
two  sections. 

2.2.2.  Plane  Channel  Flow  Stability  Equations 

The  stability  equations  for  plane  channel  flow  are  obtained  from  the  Cartesian 
Navier-Stokes  equations  Eq.  (2.19)  with  base  flow  given  by  Eq.  (2.20).  In  this  study,  the 
perturbation  is  assumed  to  be  two-dimensional.  To  investigate  the  streamwise  instability 
of  the  base  flow  to  two-dimensional  disturbances,  a  perturbation  of  the  form 
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u 


/ 


u  -  u‘e. 


(2.31) 


is  assumed,  where  Cartesian  components  of  the  velocity  mode  vector  are  used.  The  plane 
channel  linear  stability  equations  are  obtained  by  substituting  Eqs.  (2.31)  and  (2.30)  into 
the  Cartesian  Navier-Stokes  equations,  Eq.  (2.19),  and  assuming  the  base  flow  to  be  the 
analytical  solution  Eq.  (2.20).  The  following  system  of  linear  ordinary  differential 
equations  is  obtained: 


X  ^  . 

u  -  - u  -  lap  =  0 

dx^ 


—{D^-a?)+i{a>-aUd 
Re  * 


-  Dp  =  0 


ittM*  +  DCi^  =  0 


(2.32) 


where  D  represents  the  differential  operator  d  I  dx^ .  The  three  equations  can  be 
combined  into  a  single  fourth-order  ordinary  differential  equation  by  eliminating  the 
streamwise  velocity  component  and  the  pressure.  The  resulting  equation  is  called  the  Orr- 
Sommerfeld  equation: 


[{D^-a^)+  iRe (to  -  a  U^) )(D  ^-a^)+iRea 


d^U, 


(dx  ^)^ 


=  0 


or 


D*+i^e((ji  -  aU^)-2a^^D^+a^{a^ -iRe((x>  -  a  f/^)j  +i/?ca 


d^U. 


3^2 


{dx^) 


=  0 


(2.33) 
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2.2.3.  Curved  Channel  Flow  Stability  Equations 

The  stability  equations  for  curved  channel  flow  are  obtained  in  a  similar  manner. 
The  polar  Navier-Stokes  equations  Eqs.  (2.25)  and  (2.26)  are  used  with  the  base  flow 
given  by  Eq.  (2.27).  To  investigate  the  instability  of  the  base  flow  to  streamwise  two- 
dimensional  disturbances,  a  perturbation  of  the  form: 

u'  =  +  m '•(r)e‘<“®-“'>e^  (2.34) 


is  assumed  for  each  Fourier  mode.  The  curved  channel  linear  stability  equations  are: 


— (r2cg-l)+ir(r(o-aC/J 

+ 

1  .  odUu 

— 2ia  -r^ - ~fU, 

Re 

Re  dr  _ 

-—Tm  +  2rU. 

M®  + 

— (r^Sf  -  1)  +ir(rci)  -  cuf/J 

Re  \ 

Re  _ 

m'’  -  irap  =  0 
W-  r^Dp  =  0 


(2.35) 


iatZ”  +  D{ru  0  =  0 


where  +  —D  -  — 


where  D  represents  the  differential  operator  dtdr.  Again,  the  three  equations  can  be 
combined  into  a  single  fourth-order  ordinary  differential  equation  by  eliminating  the 
streamwise  velocity  component  and  the  pressure.  The  resulting  equation  is  the  curved 
channel  Orr-Sommerfeld  equation:: 
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a^D 


+  ar 


W  =  0 


where 


«3  = 

^2  =  5r^  -2a}r^  +  \Reijir‘^  -  iRear^Uf^ 
flj  =  -r-2a}r  +  3iRe(x)r^  -SiRear^Ui^ 


Uq  =  1  -2a^  +  a^  +iRe(j)r^  -iRea^oir^  -2iRearUi^  + 

-  ,dU,  d^U, 

\Rea^rU.+  iRear  ^ — -  +  [Rear  ^ - - 

dr  dr^ 


(2.36) 


See  Singer,  Erlebacher  and  Zang  (1992)  for  a  similar  derivation.  However,  in  their  study, 
the  nonlinear  evolution  equations  are  derived,  of  which  the  linear  disturbance  equation  is  a 
subset. 

2.3.  Summary 

In  §2.1,  the  3-D  Navier-Stokes  equations  were  formulated  for  the  geometry  of 
interest,  that  is  a  general  2-D  channel.  At  this  point,  the  equations  and  all  boundary 
conditions  are  prescribed  with  the  exception  of  the  initial  condition.  The  initial  condition 
will  be  obtained  using  the  linear  stability  theory  of  §2.2.  The  numerical  methodology  for 
employing  this  initial  condition  and  for  temporal  numerical  simulation  of  the  Navier- 
Stokes  equations  on  a  discrete  set  of  points  is  given  in  the  next  chapter. 


Section  3.  Numerical  Approach 


The  numerical  simulation  of  stability  and  transition  in  channel  flows  is 
accomplished  by  solving  the  continuum  governing  equations  described  in  Section  2  using 
numerical  methods.  Although  any  valid  numerical  method  for  solving  time-dependent 
partial  differential  equations  can  be  used,  spectral  methods  are  most  commonly  used  for 
solving  spatially  periodic  stability  problems.  The  standard  text  on  spectral  methods  is  that 
of  Canuto  et  al.  (1988)  and  is  the  source  of  the  information  in  this  section  relating  to 
spectral  methods.  A  detailed  example  of  the  use  of  spectral  methods  for  solving  the 
Navier-Stokes  equations  is  given  in  Zang  and  Hussaini  (1985b).  The  spectral 
discretization  used  in  periodic  stability  problems  employs  Fourier  decomposition  in  the 
periodic  directions,  thus  providing  solutions  which  are  periodic  to  the  Mh  derivative, 
where  //  -i-l  is  the  number  of  Fourier  modes  used.  The  Fourier  method  is  also  well- 
matched  to  the  temporal  simulation  problem  since  the  solutions  obtained  from  linear 
stability  theory  which  are  used  to  create  the  initial  condition  for  numerical  simulation,  are 
Fourier  modes  in  the  streamwise  and/or  spanwise  directions.  Errors  in  the  solution 
obtained  with  a  Fourier  method  can  always  be  expressed  in  terms  of  the  omitted  high 
frequency  modes  I  n  I  >  N/2  .  Therefore  low  mode-number  waves  can  be  resolved 
exactly. 

In  the  non-periodic  (wall-normal,  or  ^^)  direction,  a  Fourier  method  cannot  be 
used.  Instead,  the  solution  is  represented  in  terms  of  orthogonal  Chebyshev  polynomials, 
and  the  grid  points  are  selected  using  Gaussian  quadrature.  Though  this  method  is  more 
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computationally  intensive  per  grid  point  (or  Chebyshev  mode)  then,  for  example,  the  finite 
difference  method,  the  greater  accuracy  of  the  solution  allows  the  use  of  fewer  points  than 
other  methods.  The  basic  development  of  the  Fourier  and  Chebyshev  spectral  methods 
are  described  in  §3.1.  The  methods  are  then  applied  to  the  governing  equations  in  §3.2. 

The  discretization  in  time  is  performed  with  a  standard  second-order  finite 
difference  scheme.  The  development  of  the  temporal  discretization  and  the  iterative 
matrix  method  used  to  solve  the  discrete  system  of  governing  equations  at  each  time  step 
are  described  in  §3.3.  The  remaining  two  sections  of  this  chapter,  §§3.4  and  3.5, 
respectively,  describe  the  generation  of  the  initial  flow  solution  using  linear  stability  theory 
and  the  extraction  of  mode  information  from  the  computed  solution  for  comparison  with 
stability  theory. 

3.1.  Spectral  Methods 

The  problem  of  numerical  approximation  to  a  continuum  system  of  differential 
equations  can  be  approached  from  many  different  perspectives.  One  consistent 
methodology  is  to  formulate  a  discrete  basis  to  represent  the  spatial  variation  of  the 
solution.  In  spectral  methods,  the  discrete  basis  used  to  represent  the  solution  in  one 
direction  is  obtained  from  a  singular  Sturm-Liouville  problem  with  appropriate  boundary 
conditions.  The  Fourier  system  is  one  such  basis,  and  can  be  used  to  expand  periodic 
functions.  Other  bases  are  Chebyshev  polynomials,  Legendre  polynomials,  and  other 
Jacobi  polynomials.  These  bases  can  be  used  to  represent  non-periodic  functions.  The 
key  feature  of  spectral  methods  is  that  discretization  is  global  rather  than  local.  The 
derivatives  of  the  solution  at  a  given  point  are  determined  from  the  solution  values  at 
every  point.  Therefore,  the  order  of  accuracy  of  the  solution  increases  as  points  are 
added.  Thus  in  the  limit  as  an  infinite  number  of  discrete  points  are  used,  the  order  of 
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accuracy  becomes  infinite.  This  is  referred  to  as  spectral  accuracy  and  is  in  direct  contrast 
to  a  local  discretization  method  such  as  finite  difference,  in  which  as  the  number  of  points 
becomes  infinite,  the  order  of  accuracy  of  the  scheme  remains  fixed  at,  for  example, 
second  order. 

The  formal  arguments  as  to  the  accuracy  of  spectral  discretizations  are  more 
complex  than  those  of  finite  difference  schemes.  Refer  to  Canuto  et  a/.(1988)  for  a 
thorough  analysis  of  spectral  methods.  In  this  study,  several  numerical  experiments  are 
performed  to  demonstrate  the  accuracy  of  the  spectral  scheme.  These  are  found  in  §4.1. 

This  section  shows  an  example  of  the  use  of  spectral  methods  for  solving  a  one¬ 
dimensional  heat  equation.  Then,  the  Fourier  and  Chebyshev  discretizations  are  described 
in  detail.  The  notation  used  in  this  section  is  tailored  to  its  use  in  subsequent  sections.  In 
particular  note  that  the  Fourier  and  Chebyshev  collocation  methods  determine  the  location 
of  the  discrete  grid  points  and  give  the  derivative  operator  as  a  matrix.  This  will  be  used 
extensively  in  the  discretization  shown  in  §3.2. 

3.1.1.  Model  Problem 

To  demonstrate  the  use  of  an  orthogonal  basis  to  solve  a  partial  differential 
equation,  consider  the  1-D  heat  equation: 

(3.1) 

with  appropriate  boundary  conditions.  The  dependent  variable /(jc,r)  can  first  be  expanded 
in  terms  of  any  complete  basis  (|)„(x)  with  n  =  0,1,2.. .,<»  or  with  n  =  ,  depending 

on  the  basis.  Assuming  the  latter  is  used. 
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oo 

=  X)/„(0<i>„w 

n=  -o® 


(3.2) 


The  spatial  derivatives  of  the  basis  functions  are  needed  to  determine  the  spatial 
derivatives  of  /  .  They  must  be  expressed  in  terms  of  the  original  basis.  Once  this  is 
known,  the  next  step  is  to  expand  the  spatial  derivative  of/ in  terms  of  the  orthogonal 
basis.  In  this  case,  for  example, 

—  (|>„(^)  =  E  (3.3) 

dx  m=  -oo 


Once  this  is  known,  the  next  step  is  to  expand  the  spatial  derivative  of/ in  terms  of  the 
orthogonal  basis. 


— /(^.O  =  E/„(0— -  =  E/„(0  E 

dx  n=-oo  dx 


n  =  -oo  m=-«> 


^  E  where  ^  Smpfjt) 

m  =  -oo 


p  =  -<^ 


(3.4) 


A  spatially  discrete  system  is  then  easily  obtained  by  truncating  the  expansion  to  include 
only  A^+1  terms  n  =  -  M2,...,  M2  .  Substituting  the  discrete  expansion  of  /and  its  second 
derivative  into  Eq.  (3.1)  yields  a  system  of  M+1  ordinary  differential  equations  in  time.  In 
the  example. 


(3.5) 
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Once  the  relationship  between  the  basis  functions  and  their  derivatives  is  determined,  that 
is,  the  constants  ,  most  of  which  are  zero,  the  system  can  be  easily  solved. 

The  goal  of  the  spectral  discretization  is  to  transform  the  original  differential 
equation  into  +1  ordinary  differential  equations  in  time.  The  boundary  conditions  are 
implemented  in  various  ways  depending  on  the  type  of  scheme  used.  The  crux  of  the 
discretization  is  the  determination  of  the  spatial  derivatives  of  the  basis  functions.  Then,  a 
self-consistent  set  of  coupled  ordinary  differential  equations  for  the  amplitude  of  each 
basis  function  is  obtained. 

Two  approaches  to  setting  up  this  system  of  ODEs  are  used  in  this  study.  The 
Fourier  Galerkin  method,  described  in  the  next  section,  uses  a  Fourier  basis  to  represent 
the  solution  to  the  governing  equations.  The  equations  themselves  are  transformed  into 
Fourier  space.  The  second  and  third  are  both  called  collocation  methods.  These  methods 
use  the  basis  function  formulation  to  derive  a  method  of  differentiating  the  solution, 
without  the  need  to  transform  the  solution  into  Fourier  or  Chebyshev  space.  The 
derivative  is  obtained  as  a  matrix  multiplication  operator. 

3.1.2.  Fourier  Galerkin  Method 

The  Fourier  Galerkin  spectral  method  uses  a  trigonometric  polynomial  basis,  or 
Fourier  basis,  to  represent  the  solution  to  the  governing  equations.  The  application  of  the 
Fourier  basis  to  the  above  model  problem  is  straight-forward.  A  domain  ofjc  e  [0,2iT]is 


assumed. 
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<1>„W  = 

00 

«  =  -«> 

-^f(x,t)  =  E  -«Y„(0e‘""  =  E/„.;^(0e‘ 

OX  n  =  -«»  «  =  -<» 

where  /„^J0  =  -n^f^(t) 


(3.6) 


and  the  resulting  system  of  ordinary  differential  equations  would  be,  after  truncation: 

,  n=-N/2,...,NI2  (3.7) 

The  Fourier  basis  has  the  unique  property  that  derivatives  of  a  basis  function  are 
proportional,  with  a  complex  constant,  to  the  basis  function.  Thus  the  ODEs  that  are 
obtained  from  discretization  are  decoupled.  Notice,  however,  that  this  occurs  only  for 
linear  differential  equations  with  coefficients  that  are  constant  in  the  direction  of 
discretization!  Nonlinear  terms  and  non-constant  coefficient  terms  introduce  coupling  in 
the  ODEs.  However,  these  terms  will  not  be  treated  using  a  Fourier  Galerkin  method  in 
this  study,  but  will  instead  employ  the  Fourier  collocation  method  of  the  next  section. 
3.1.3.  Fourier  Collocation  Method 

The  Fourier  collocation  method,  and  the  Chebyshev  collocation  method  of  the 
following  section  use  orthogonal  bases  to  satisfy  the  governing  equations  as  did  the 
Fourier  Galerkin  method.  However,  the  Fourier  Galerkin  method,  by  solving  the 
governing  equations  in  transform  space,  assumes  that  the  solution  satisfies  the  transformed 
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governing  equations  everywhere  in  the  domain,  and  that  the  solution  will  have  an  error  of 
the  order  of  the  lowest  frequency  mode  neglected  in  the  approximation.  However, 
collocation  methods  assume  that  the  basis  functions  hold  only  at  certain  points,  called  the 
collocation  points,  which  are  given  by  a  quadrature  rule.  Thus  the  solution  is  represented 
as  iV +1  discrete  values  rather  than  N  +1  mode  amplitudes.  The  basis  is  used  to  determine 
the  relationship  between  the  solution  and  its  derivatives. 

In  the  Fourier  collocation  method,  again  trigonometric  polynomials  are  used  as  a 
basis.  The  solution  is  determined  at  N+l  collocation  points,  where  is  A^is  even.  The 
Gauss-Lobatto  collocation  points  are  used: 


271/ 

N+l 


i=0,...,N 


(3.8) 


where  the  dependent  variables  are  periodic  in  ^  with  period  27T  and  are  represented  on  the 
computational  domain  ^  G  [0,27t].  The  values  of  the  dependent  variables  at  the  collocation 
points  are  given  by  the  expansion,  (usually  called  the  inverse  Fourier  transform): 

N/2 

=  E  (3.9) 

k=-N/2 

and  the  expansion  coefficients  are  given  by  the  Fourier  transform: 

A  =  ^  E  /(^,)  (3.10) 

W+l  /=o 


Differentiation  of  the  dependent  variables  is  accomplished  through  the  expansion: 


50 


M. 


NI2 

E  ii/.e'*'- 

k  =  -m 


(3.11) 


The  transformations  and  differentiation  can  be  expressed  strictly  as  matrix  operations. 
First,  the  dependent  variable  and  its  transform  are  expressed  in  vector  notation: 

F  =  vector  {/(^;)  |/  =  0,Ar} 

^F.  vector  {^«,)|i=0,Wl  (3., 2) 

F  H  vector  {f^\k=  -N/2 ,N/2} 


The  Fourier  transform  and  inverse  transform  are  then  expressed  as  matrix  multiplies: 

F  =  Cp^F 

where  {Cp^]j^  =  e‘*^^  , 

eind  differentiation  is  also  represented  as  a  matrix  multiply: 


F  =  CpF 


(3.13) 


N+1 


-ik^j 


—F  =  D^F 
dl  ^ 

Dp  =  Cf^'diag{i*}C^ 


(3.14) 


The  matrix  derivative  Dp  can  also  be  evaluated  in  closed  form  using: 


1 

iV+1 


m 


E 

k=-m 


2^(/-y)7t^ 

AT+l  j 


(3.15) 
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The  model  problem,  Eq.  (3.1) ,  can  be  formulated  using  the  Fourier  collocation  method. 
Using  the  notation  of  Eq.  (3.12),  the  heat  equation  can  be  expressed  as: 


dt 


=  -dIf 


(3.16) 


Also,  a  nonlinear  problem  may  also  be  formulated  using  the  Fourier  collocation 
method.  For  example  the  nonlinear  Burger's  equation. 


^/(x,r)  =  -f{x,t) 
dt 


df(x,t) 

dx 


(3.17) 


with  periodic  boundary  conditions  can  be  discretized  spatially  as: 


dF 

dt 


-F  *  (D^F) 


(3.18) 


where  the  multiplication  in  the  nonlinear  term  is  a  pointwise  vector  multiply,  denoted  by 
the  symbol  *.  In  the  remainder  of  this  document,  pointwise  multiplication,  rather  than  an 
inner  or  outer  product,  is  assumed  when  two  vectors  are  multiplied  and  the  symbol  *  is 
omitted. 

To  review,  a  single  matrix  operator  D^-  is  derived  which  represents  the 
combination  of  transforming  the  solution  into  Fourier  space,  taking  the  derivative  by 
scaling  each  expansion  coefficient,  and  transforming  the  solution  back  into  physical  space. 
Note  that  the  periodic  boundary  condition  on  the  solution  is  automatically  satisfied  by  the 
choice  of  basis  functions  and  that  no  additional  boundary  conditions  are  needed.  See 
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Canute  et  a/.(1988)  for  more  details  on  the  Fourier  collocation  method,  including 
estimates  of  accuracy. 

3.1.4.  Chebyshev  Collocation  Method 

In  the  Chebyshev  collocation  method,  a  Chebyshev  polynomial  basis  is  used  to 
represent  the  solution.  Chebyshev  polynomials  are  given  by: 

=  cos(mcos-'(^.))  (3.19) 

Trigonometric  identities  give  the  following  recursion  relations  for  Chebyshev  polynomials: 

r„®=i 

J'l®'?  (3.20) 

r..,(0  =  25r„(5)-7;.,(0  ,  mil 

The  derivative  of  a  Chebyshev  polynomial  is  found  in  terms  of  other  Chebyshev 
polynomials  using  the  recursion  properties.  The  recursion  mle  for  the  derivative  is: 

27'.®  =  ^t.i®  -  ^C,®  «’il  (3.21) 

m+1  m-1 


To  discretize  a  differential  equation  using  the  Chebyshev  collocation  method,  first,  the 
collocation  points  are  defined  from  a  Chebyshev  quadrature  rule: 


=  -cos 


^  — 
'n 


(3.22) 


Then,  the  solution  at  the  collocation  points  is  expanded  in  terms  of  the  Chebyshev  basis: 

/®)=  E/.r.®) 

m=-N 


(3.23) 
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Again,  the  solution,  its  transform,  and  its  derivative  are  represented  as  solution  vectors  of 
length  AT +1: 


F  =  vector  {/(^,.)  |i  =  0,A} 
—F  =  vector  { — (^.)  |/  =  0,A} 

F  =  vector  {/  |m  =  0,A} 


(3.24) 


The  Chebyshev  transform  and  inverse  Chebyshev  transform  are  then  given  as  matrix 
multiplication  operators; 


F  =  CJ'f 


F  =  C,F 


{CJ  =  cos 


N) 


{Ccb.  = 


Nc.c, 


COS 


I 


(3.25) 


c  =  < 

m 


2  m  =  0,N 
1  l^m<A-l 


Differentiation  is  again  accomplished  through  a  matrix  multiplication: 


—F  =  D^F 
Dc  ~  ^ 


(3.26) 


The  R  operator  above  refers  to  the  application  of  the  recursion  relationship  of  Eq.  (3.21), 
which  can  be  expressed  as  a  matrix  multiply.  As  in  the  Fourier  collocation  method,  a 
closed  form  exists  for  the  Chebyshev  derivative  operator: 
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{^r}  = 

L  *  mn 


(-1) 


m+n 


C  ^  ^ 


-L 


m  *  n 


1  <  m=n  <  N 


2(1-0 

2N^+l 

6 

2N^+l 


m  =  n  =  1 


m  =  n  =  N 


(3.27) 


The  Chebyshev  collocation  method  can  be  applied  to  the  above  model  problems. 
However,  the  boundary  conditions  cannot  be  periodic.  Rather,  say  the  end  values  are 
fixed,  and  the  domain  is  given  as  xe  [-1,1] .  The  first  problem  is  discretized  as; 


dF 

dt 


=  -dIf 


(3.28) 


The  boundary  conditions  are  imposed  directly  by  removing  the  two  equations  at  the 
boundaries  and  replacing  them  with  conditions  that  the  statement  of  fixed  boundary 


conditions: 


^0  -  /(-I) 
^v=/(l) 


(3.29) 
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3.1.5.  Chebyshev  Collocation  Method  on  a  Staggered  Grid 

In  this  study  a  staggered  grid  is  needed  in  the  direction  in  which  Chebyshev 
polynomials  are  used.  Specifically,  the  pressure  solution  and  continuity  equation  are 
located  at  points  staggered  from  the  velocity  solution  and  momentum  equation.  This  can 
be  easily  accomplished  within  the  collocation  framework.  The  "half-cell"  grid  points, 
called  the  half  points  for  lack  of  a  better  name,  are  prescribed  by  the  quadrature  rule: 


(3.30) 


The  Chebyshev  transform  and  inverse  transform  on  the  half  points  are  modified  slightly: 


.1 

(C,  =  cos 


2  T^(j+-)k 

iCch,  =  cos  — ^ 

NCjC^  \  N 


(3.31) 


for  0^7^(A^-1)  and  0<k<N 


The  operators  on  the  half  points  are  denoted  with  a  “o”.  The  grid  points  are 
denoted  with  “+”.  Interpolation  between  the  grid  points  and  the  half  points  is 
accomplished  by  Chebyshev  interpolation.  To  interpolate,  for  example,  the  pressure  P 
from  the  half  points  to  the  grid  points,  the  pressure  is  transformed  from  the  half  point 
physical  values  into  Chebyshev  space,  and  then  from  Chebyshev  space  onto  the  regular 
grid  points.  This  can  be  represented  as  a  matrix  operation: 

p "  =  a;  P  " 
where 


(3.32) 
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where  A  "  represents  interpolation  from  the  grid  points  to  the  half  points.  The  inverse  of 
this  matrix  is: 

a;  =  (a;)-'  =  {c^r^Cc  0.33) 

and  represents  the  reverse  interpolation. 


3.1.6.  Chebyshev  Quadrature 

To  determine  the  kinetic  energy  of  the  disturbance  in  the  simulation,  an  integral  is 
needed  across  the  Chebyshev  grid.  This  can  be  accomplished  using  Chebyshev 
quadrature.  The  weights  for  Chebyshev  quadrature  are  computed  by  integrating  the 
recursion  rule  Eq.  (3.21).  The  integration  is  given  by: 


-1 


w. 


^  4c 


0  N(l-m^)  N 

1 


cos(-^^)  0<i<N 


(3.34) 


i=0,A^ 


3.1.7.  Summary 

The  use  of  basis  functions  in  all  three  spatial  directions  allows  discrete  derivative 
operators  to  be  derived,  assuming  that  the  solution  is  represented  on  certain  collocation 
points.  In  the  next  section,  these  operators  and  collocation  points  are  used  to  discretize 
the  continuum  Navier-Stokes  equations  formulated  for  general  2-D  channels. 
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3.2.  Spatial  Discretization  of  Governing  Equations 

The  general  approach  for  the  spatial  discretization  of  the  Navier-Stokes  equations 
for  a  generalized  2-D  channel  geometry  is  to  use  1)  a  Fourier  Galerkin  method  in  the 
spanwise  direction,  excluding  the  nonlinear  terms,  2)  a  Fourier  collocation  method  in  the 
streamwise  direction  and  to  evaluate  the  nonlinear  terms  in  the  spanwise  direction,  and  3) 
a  Chebyshev  collocation  method  in  the  wall-normal  direction.  The  mixed  use  of  Fourier 
Galerkin  and  Fourier  collocation  in  the  spanwise  direction  is  consistent  with  the  literature, 
but  note  that  in  some  of  the  literature  this  type  of  method  would  be  called  a  "pseudo- 
spectral  method"  or  a  would  retain  the  name  "collocation  method."  In  both  cases,  the 
implication  is  that  the  method  is  not  truly  Galerkin  if  some  terms  are  evaluated  using  a 
collocation  method. 

The  approach  to  the  spatial  discretization  in  this  study  is  similar  to  that  used  by 
others  to  study  plane  and  curved  channel  flows.  See  the  review  of  Kleiser  and  Zang 
(1991)  and  the  simulations  of  Orszag  and  Patera  (1983),  Finlay,  Keller  and  Ferziger 
(1988),  Streett  and  Hussaini  (1991)  and  Zang  and  Hussaini  (1987).  However,  this  study 
differs  notably  from  those  above  in  that  a  general  curvilinear  coordinate  system  is  used  in 
the  plane  of  the  streamwise  and  wall-normal  directions.  Two  examples  of  studies  of  2-D 
channel  flows  incorporating  general  curvilinear  coordinates  are  Kamiadakis,  Mikic  and 
Patera  (1988)  and  Guzman  and  Amon  (1994),  both  of  which  used  the  spectral  element 
numerical  method.  This  method  is  quite  different  from  the  spectral  method  used  here. 

The  spectral  element  method  is  closely  related  to  the  finite  element  method  and  uses  local 
basis  functions  to  represent  the  solution  quantities  on  small  geometric  elements.  The 
present  method  uses  global  basis  functions  as  described  in  the  previous  chapter. 
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The  discrete  form  of  the  Navier-Stokes  equations  will  be  derived  first  by 
transforming  the  equations  into  Fourier  space.  The  nonlinear  terms  are  included  as  a 
source  term  which  will  be  computed  using  the  Fourier  collocation  method.  Second,  the 
solution  will  be  discretized  on  the  Fourier-Chebyshev  collocation  grid,  with  pressure  and 
continuity  staggered.  The  continuum  derivative  operators  will  be  replaced  with  discrete 
spectral  operators.  In  addition,  the  metrics  in  the  equations  will  be  replaced  with  discrete 
vectors  containing  values  of  the  metrics  at  each  point. 

3.2.1.  Discretization  in  the  Spanwise  Direction 

The  contravariant  component  momentum  equations,  Eq.  (2. 15),  are  Fourier 
transformed  in  the  spanwise  direction: 
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for  ^  =  -  N 2I2  ,  N2I2,  and  m  -  1,2,3.  The  nonlinear  convective  term  is  given  by  F^,  , 
which  represents  the  k  mode  of  the  transformed  nonlinear  term.  The  discretization  of  the 
nonlinear  term  in  physical  space  is  given  in  Eq.  (3.42). 

The  transformed  continuity  equation  in  terms  of  contravariant  velocities  is: 


fgdl 


+  ikut 


=  0 


(3.36) 


3.2.2.  Discretization  in  Streamwise  and  Wall-Normal  Directions 

The  spatial  differentiation  in  the  streamwise  and  wall-normal  directions  is 
accomplished  through  the  Fourier  and  Chebyshev  collocation  methods,  respectively, 
described  above.  The  discrete  solution  vector  is  represented  as: 


U. 


=  vector {M*'"(^f,^J)  |  ,y  =  0,A^3} 

P"  =  vector{p™(^;,^'  )  I  1  =  0, IVi  ,y  =  0,lV3-l} 

2 


(3.37) 


for  ^  =  -  N2/2  ,  N2/2.  The  metric  terms  are  written  as  diagonal  matrices  so  that  matrix 
multiplication  can  be  used  to  represent  the  products  in  Eqs.  (3.35)  and  (3.36).  These  are 
written  as: 


(3.38) 
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To  obtain  the  spatially  discrete  governing  equations,  derivatives  in  the  streamwise 
and  wall-normal  directions  are  replaced  by  their  matrix  equivalents.  In  addition,  a  matrix 
collocation  derivative  is  required  in  the  spanwise  direction  for  the  evaluation  of  the 
nonlinear  term.  These  matrix  operators  are: 


0 

D.  =  —  =  Fourier  Collocation  in 
' 

0 

£>2  =  — -  =  Fourier  Collocation  in  ^  ,  for  use  in  non-linear  term 

se 

0 

Dj  =  —  =  Chebyshev  Collocation  in 


(3.39) 


Additional  multiplications  by  Chebyshev  interpolation  matrices  are  required  to 
account  for  the  staggered  pressure  and  continuity  equation.  The  spatially  discrete 
momentum  equations  are: 


dt 


Re  * 


y(£,D,+£,D3)(?;  .  r;Hfil .  s;ui 


Re 


*  Fj't  =  -K  -  *  Q- 

where 

E.  =  G-\D^GG‘^  +  D^GG‘^) 

F.  =  G^'D^  +G‘^D^ 

H,  =  +  F,. 


and  P,  =\P, 


(3.40) 
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ior  k  =  -  N 2I2,  N 212  ,  and  m=  1,2,3,  where  the  nonlinear  convective  terms  are  calculated  in 
physical  (non-transformed)  space  and  then  transformed  into  Fourier  space.  The  solution 
vector  representations  of  the  velocity  and  nonlinear  term  are: 

f/'"  =  vector I  i=0,N,  ,  ,  j=0,N,}  = 

f/"”  =  vector{M7(^;,^;)  I  k=-N^/2,N^2  ,  i=0,N^  J=0,N,} 

9-'^  =  vector{  I  i=0,N^  ,  n=0,N2  ,  ;=0,A^3}  =  (3.41) 

=  vector I  k=-N^/2,N^/2  ,  /=0,iV,  ,  j=0,N^} 

where  f  =  =  V-(««) 


Using  the  above  notation,  the  nonlinear  term  is  written: 

=  JiU't/"'  +  + 

where 

7,  =  G'DjG 

72  =  G'D^G 

73  =  G'^DjG 


(3.42) 


The  Fourier  transforms  are  accomplished  through  matrix  multiplications.  The  steps  for 
computing  the  nonlinear  term  are: 

1)  G'”  =  c;'  u"' 

2)  Compute  (3.43) 

3)  ^  =  Cp 


Finally,  the  continuity  equation,  interpolated  onto  the  staggered  grid,  is: 

A-(d,GU;  .  ikU^  .  =  0 


(3.44) 
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3.2.3.  Final  Spatially  Discrete  Matrix  System 

The  above  system  of  equations  must  be  solved  in  a  coupled  manner  in  order  to 
determine  the  discrete  solution  at  all  times.  The  coupled  matrix  system,  still  continuous  in 
time,  can  be  concisely  expressed  as: 


for  ^  =  -  N2I2,  N2/2.  The  terms  L  ,G ,  and  D  refer  to  the  linear  Laplacian,  gradient,  and 
divergence  coefficients,  respectively,  obtained  from  placing  Eqs.  (3.40)  and  (3.44)  in 
matrix  form.  The  A^^+l  systems  are  coupled  through  the  term^R^^,  which  represents  the 
three  terms  on  the  right  hand  side  of  Eq.  (3.40) 

Boundary  conditions  are  implemented  directly  by  omitting  equations  from  the 
above  matrix  system  and  replacing  them  with  implicit  statements  of  the  boundary 
conditions.  The  periodic  condition  in  the  streamwise  and  spanwise  directions  is 
automatically  satisfied  by  the  Fourier  discretization.  At  the  walls,  the  velocity  is  specified. 
While  no  boundary  conditions  on  the  pressure  are  required,  it  is  necessary  to  provide  one 
pressure  condition.  Since  the  pressure  appears  in  the  governing  equations  only  in  terms  of 
its  gradient,  the  pressure  is  determined  only  up  to  an  arbitrary  constant.  This  is  manifested 
in  the  current  formulation  as  a  single  zero  eigenvalue  in  the  linear  coefficient  matrix,  A  , 
due  to  a  redundant  equation  in  the  continuity  equation.  Therefore,  one  continuity 
equation,  at  an  arbitrary  location,  is  removed  and  replaced  by  a  condition  which  sets  the 
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arbitrary  value  of  pressure.  This  is  achieved  by  requiring  that  the  mean  value  of  pressure 
be  zero.  After  taking  the  Fourier  transform,  this  condition  becomes: 

E  E  w-  =  0  (3.46) 

1=0  i=0 

for  each  value  of  k,  where  Wj  represent  the  Chebyshev  quadrature  weights  of  Eq.  (3.34). 

For  a  given  discrete  initial  condition,  the  above  system  of  discrete  equations,  with 
prescribed  wall  velocities  and  the  mean  pressure  condition,  is  a  well-posed  set  of  coupled 
ordinary  differential  equations  in  time.  It  does,  however,  have  the  problem  of  being  highly 
ill-conditioned,  as  will  be  discussed  in  Section  4. 

3.3.  Temporal  Discretization  and  Solution  Method 

The  solution  of  a  set  of  linear  first-order  coupled  ordinary  differential  equations 
(ODEs)  is  a  commonplace  problem  in  numerical  analysis.  The  system  is  linear  in  the  time 
derivative,  but  has  a  coupled  nonlinear  function  as  its  source  term.  This  type  of  problem 
can  be  solved  using  any  of  a  number  of  standard  procedures.  Explicit  or  implicit  methods 
of  various  orders  can  be  used.  For  this  study,  a  second-order  implicit  trapezoidal 
integration  was  selected.  Iteration  is  required  at  each  time  step  to  update  the  nonlinear 
terms.  The  scheme  could  easily  be  replaced  with  a  higher-order  implicit  method. 

Though  the  discretized  system  is  a  system  of  ODEs,  the  partial  differential 
equation  (PDE)  nature  of  the  original  equations  has  a  severe  effect  on  the  nature  of  the 
coupled  discrete  system.  In  particular,  the  continuity  equation  acts  as  an  instantaneous 
compatibility  equation  between  the  velocity  components  which  must  be  satisfied  at  every 
point.  The  approach  used  here  is  that  while  the  discretized  Navier-Stokes  equations  can 
be  solved  using  any  of  a  number  of  standard  procedures  for  linear  first-order  coupled 
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ODEs,  the  tight  coupling  enforced  by  the  continuity  equation  must  be  carefully  accounted 
for.  For  this  reason,  the  continuity  equation  must  be  completely  implicit  at  each  time  step 
and  for  each  iteration.  Thus  the  tight  coupling  between  the  solution  components  will  be 
maintained  even  as  the  nonlinear  terms  are  iteratively  updated.  Numerical  experiments  in 
which  continuity  was  solved  iteratively  all  showed  rapid  divergence. 

The  application  of  the  trapezoidal  integration  scheme  to  the  spatially-discrete 
Navier-Stokes  equations  is  shown  in  the  following  section,  while  the  iterative  procedure 
for  advancing  the  solution  from  one  time  level  to  the  next  is  described  in  §3.3.2. 

3.3.1.  Temporal  Discretization 

The  integration  of  the  spatially  discrete  system  of  equations,  Eq.  (3.45),  can  be 
performed  using  any  numerical  method  for  solving  initial- value  ordinary  differential 
equations.  In  this  study  a  second-order  central  difference  discretization  is  applied.  Time 
is  discretized  with  a  constant  time  step,  t  =  nAt.  The  entire  equation  is  evaluated  at  the 
half  time  level  n-i-1/2.  This  is  achieved  by  averaging  the  solution  at  the  n  and  n+1  time 
levels.  The  discretized  form  of  Eq.  (3.45)  is: 


Vrn+l 

+  u, 

1 

2 

/  N 

£  1 

►  + 

H  2 

At 

0 

^/i+l 

k 

0 

2 

This  is  re-written  in  terms  of  the  change  in  the  solution  vector  at  each  time  step.  The  final 
discrete  version  of  the  governing  equations  is  given  in  two  forms: 
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f  ^  ^ 

v:  *  |a(7, 

*  +  A 

►  =  ^ 

2  /  . 

[  ®  J 

0 

where 


Au,-ir;'-u: 

Au  = 


(3.48) 


^  =  At 

-A  ^ 

^  fl 

K 

^  ^ 

0 

where 


I  +  ALl 
2  2 

-D  0 

2 


(3.49) 


In  the  second  form,  the  linear  matrix  B  is  multiplied  by  the  vector  of  unknowns.  The  first 
term  on  right  hand  side  is  known  from  the  previous  time  level,  while  the  second  term  is 
dependent  on  the  unknown  AU.  Thus  the  system  of  equations  is  of  the  form  LX  =  N(X), 
where  X  is  the  vector  of  unknowns,  L  is  a  linear  coefficient  matrix,  and  N(X)  is  a  nonlinear 
function  of  X  An  iterative  method  is  needed  in  order  to  solve  this  system. 


66 


3.3.2.  Iterative  Solution  Method 

The  matrix  system  Eq.  (3.49)  is  solved  at  each  time  step  using  an  iterative  solution 
technique.  At  each  time  step,  first  the  change  in  the  solution  vector  is  guessed  as 
Al/j  =  AI7°  ,  AP^  =  AP°.  The  above  system  is  then  used  to  determine  a  correction  to 
the  guess: 


where 


At  ■  Residual 


Residual  =  -  —  B  ' 

\vl 

Rjf/"+-!-Al/)l 

1  1 

^  -  A  ' 

1  1 

►  H-  ^ 

n  2  / 

At 

<1 

_ J 

0  1 

AU, 

'au: 
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►  =  i 

>  +  i 

► 

AP- 

(3.50) 


where  B  is  some  approximation  to  the  original  B.  This  approximation  gives  a  certain 
degree  of  flexibility  in  the  iterative  method.  For  example,  if  the  choice  B  -  I  is  used,  the 
method  collapses  to  point-Jacobi.  The  iteration  is  repeated  until  the  maximum  absolute 
value  of  the  residual  is  below  some  convergence  criterion.  Note,  however,  that  the 
continuity  equation  is  unstable  to  the  above  iterative  technique,  so  a  minor  revision  is 
necessary.  The  iterative  system  is  actually  solved  by  requiring  that  all  corrections  to  the 
velocity  vector  be  divergence  free;  that  is,  every  velocity  solution  and  correction  is 
required  to  satisfy  the  continuity  equation.  With  this  modification,  the  final  iterative 
system  is  written: 
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B 


r  =  Lt  {  I  0}'  Residual 


(3.51) 


where  B  contains  some  approximation  to  the  momentum  equation  operator,  but  contains 
the  exact  continuity  equation  operator.  The  continuity  equation  has  been  removed  from 
the  right  hand  side  (by  multiplication  by  zero). 

To  summarize  the  solution  procedure,  at  each  time  step  the  change  in  the  solution 
is  guessed.  This  guess  is  selected  as  the  change  in  the  solution  at  the  previous  time  level. 
Using  the  guessed  solution,  the  residual  of  the  momentum  equation  is  computed.  The 
maximum  absolute  value  of  the  residual  at  each  iteration  is  computed  and  compared  to 
some  convergence  criterion  e.  For  the  results  in  this  study,  the  convergence  level  was  e 
=  10’*^ .  Typically,  2  to  8  iterations  were  needed,  depending  on  the  size  of  the  time  step. 
For  each  iteration,  the  update  to  the  solution  was  determined  using  the  block  Lower- 
Upper  (LU)  decomposition  method.  The  block  LU  method  is  a  version  of  direct  block 
Gaussian  elimination  in  which  the  coefficient  matrix  is  factored  into  lower  triangular  and 
upper  triangular  block  matrices,  and  the  solution  is  found  by  forward  elimination  and  back 
substitution  of  the  right  hand  side  into  the  factored  coefficient  matrix.  The  factorization 
need  be  performed  only  once  since  the  coefficient  matrix  B  does  not  depend  on  the 
solution  quantities.  The  matrix,  B  ,  which  can  be  any  suitable  approximation  to  B,  is 
selected  by  removing  terms  in  the  Laplacian  so  that  does  not  depend  on  the  spanwise 
wavenumber,  k.  The  same  B,  computed  and  factored  one  time,  is  used  for  all  of  the 
spanwise  modes.  Even  with  this  choice  of  R,  it  consumes  approximately  75%  of  the 
memory  needed  for  the  simulation. 
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3.3.3.  Comparison  with  Previous  Studies 

The  direct  approach  taken  here,  i.e.,  an  implicit  iterative  method  for  the  solution  of 
the  fully-coupled  discretized  equations  of  motion,  is  somewhat  unique.  Previous  channel 
flow  simulations  ( i.e.,  Zang  and  Hussaini,  1985, 1987,  1990)  have  used  methods  for 
decoupling  the  momentum  equation  from  the  continuity  equation  by  solving  a  Poisson 
equation  for  pressure.  They  also  do  not  solve  the  entire  equation  implicitly;  the  nonlinear 
terms  are  lagged  at  each  time  level  using  a  third-  or  fourth-order  Runge-Kutta  method. 
This  removes  the  need  for  iteration  as  is  used  in  this  study.  See  Canuto  et  al.  (1988)  for  a 
complete  discussion  of  the  solution  methods  used  by  previous  authors. 

The  advantages  of  an  implicit  iterative  method  over  an  implicit/explicit  method  are 
mostly  heuristic.  It  is  satisfying  to  know  that  the  residual  of  the  discrete  system  will  be 
forced  below  a  convergence  criterion  at  each  time  step.  This  is  not  done  in  partially 
explicit  methods.  The  numerical  stability  properties  of  implicit  methods  are  somewhat 
better  than  explicit  methods,  though  the  iteration  introduces  its  own  limits  on  the 
maximum  allowable  time  step.  Finally,  the  decoupling  of  the  momentum  and  continuity 
equations  used  in  previous  studies  requires  the  use  of  approximate  boundary  conditions  on 
pressure.  The  validity  of  various  approximate  boundary  conditions  for  pressure  has  been 
debated.  It  is  certain  that  any  approximate  boundary  condition  imposed  on  the  pressure 
will  introduce  some  error  in  the  flow  solution,  bu  tit  is  argued  that  these  errors  are  of  the 
order  of  truncation  error.  Again,  see  Canuto  et  al.  (1988)  for  details.  In  this  study,  the 
momentum  and  continuity  equations  are  solved  fully  coupled,  removing  the  need  for  such 
approximations. 

The  present  implicit  method  is  an  efficient  choice  for  the  general  curved  channel 
problem.  Though  the  memory  requirements  are  large  (approximately  64  megabytes  for  a 
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8x64x8  grid),  today’s  supercomputers  have  sufficient  memory  for  such  computations.  In 
addition,  the  complexity  of  the  general  coordinate  system  primarily  arises  in  the 
computation  of  the  Laplacian  term  in  A  (and  therefore  in  B  and  B).  This  term  is 
computed  only  once  at  the  beginning  of  the  simulation.  Finally,  the  forward  elimination 
and  back  substitution  step  required  each  iteration  is  quick.  It  uses  low-level  matrix 
operations  which  are  easily  vectorized  or  parallelized  to  take  advantage  of  the 
supercomputer  architecture.  In  summary,  the  numerical  approach  is  efficient  and 
appropriate  for  the  general  channel  problem. 

3.4.  Initial  Condition 

The  spatially  and  temporally  discrete  system  Eq.  (3.49)  is  complete  and  well- 
posed.  Periodicity  is  ensured  in  two  spatial  dimensions,  and  the  velocity  boundary 
conditions  are  implemented  directly,  while  no  boundary  conditions  are  needed  for 
pressure.  The  only  remaining  aspect  of  the  problem  is  the  initial  condition.  In  a  temporal 
simulation  of  channel  flow  instabilities,  the  initial  condition  is  of  utmost  importance.  If  an 
initial  condition  is  specified  as  simply  the  base  flow  with  no  disturbance,  the  numerical 
method  would  not  provide  enough  truncation  error  or  round-off  error  to  trigger  any 
instabilities  in  a  finite  time;  since  a  spectral  method  is  being  used,  the  errors  would  be  of 
the  order  of  machine  accuracy.  Thus  while  in  some  studies,  perturbations  which  trip 
instabilities  are  assumed  to  exist  due  to  numerical  error,  this  is  never  the  case  in 
simulations  of  transition,  in  which  the  error  is  controlled  to  be  very  small,  and  any 
disturbances  to  the  base  flow  must  be  carefully  prescribed. 

The  initial  condition  in  a  temporal  simulation  is  the  entry  point  for  information 
obtained  from  linear  stability  theory.  For  prescribed  dimensions  of  the  geometry  and  for 
prescribed  Reynolds  number,  linear  stability  theory  is  used  to  give  disturbance  waves  to 
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start  the  nonlinear  simulation.  As  was  mentioned  in  Section  1,  in  the  physical  problem 
under  consideration,  it  is  assumed  that  receptivity  mechanisms  have  tripped  the  primary 
linear  instability  mode,  so  that  it  has  grown  to  a  perceptible  amplitude.  The  actual 
amplitude  of  the  linear  wave,  however,  is  critical  to  its  evolution.  Nonlinear  effects  occur 
only  for  disturbances  of  sufficiently  large  amplitude.  It  is  conceivable  that  in  a  particular 
experiment  or  realistic  flow  situation,  that  the  disturbance  environment  will  vary.  Thus,  in 
this  study,  perturbations  of  various  initial  amplitudes  are  examined,  from  well  within  the 
linear  regime  to  very  large  amplitude  disturbances  for  which  even  weakly-nonlinear 
stability  theory  is  not  valid. 

The  initial  condition  is  comprised  of  the  base  flow  plus  a  primary  instability  wave 
at  a  given  amplitude,  plus  other  higher  modes  of  various  amplitudes.  These  higher  modes 
are  also  obtained  from  linear  theory  as  harmonics  of  the  primary  mode.  Though  the  initial 
amplitude  of  these  harmonics  proved  to  be  of  little  importance  in  this  study,  as  is  shown  in 
Section  4,  they  are  included  for  completeness  and  could  have  important  effects  in  future 
studies. 

The  remainder  of  this  section  describes  the  method  for  determining  the  initial 
condition  from  the  results  of  linear  stability  theory.  The  initial  conditions  for  plane  and 
curved  channel  flow  with  streamwise  disturbances,  used  to  obtain  the  results  of  Section  4, 
are  presented.  The  development  of  other  specific  types  of  initial  conditions,  such  as 
spanwise  modes  for  the  plane  and  curved  channel,  and  initial  conditions  for  other  channel 
shapes  will  be  reserved  for  the  next  phase  of  this  work. 

3.4.1.  Discrete  Linear  Stability  Modes  for  Plane  Channel  Flow 

To  study  the  stability  of  plane  or  curved  channel  flow,  linear  stability  theory  is  used 
to  generate  a  perturbation  to  the  base  flow.  The  base  flow  plus  perturbation  is  then 
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supplied  as  the  initial  value  of  the  solution.  As  the  solution  is  marched  in  time  using  the 
above  method,  the  evolution  of  the  perturbation  is  ascertained.  The  perturbation  is 
selected  as  the  least  stable  eigenmode  of  the  linearized  stability  equations.  The  eigenmode 
profiles  were  obtained  using  the  linear  stability  code  of  Herbert  (1988a)  using  the  fourth- 
order  Orr-Sommerfeld  equation,  Eq.  (2.33).  This  code  was  used  to  compute  the  wall- 
normal  velocity  component  mode  shape  and  to  provide  the  corresponding  eigenvalue 
(complex  temporal  growth  rate)  for  comparison  with  numerical  results.  Given  the  wall- 
normal  mode  shape,  the  streamwise  mode  shape  and  pressure  were  computed.  The 
streamwise  velocity  was  computed  directly  from  the  continuity  equation,  while  the 
pressure  mode  was  computed  from  the  streamwise  momentum  equation. 

For  plane  channel  flow,  the  mode  shapes  were  calculated  using: 


m'(.x^)  =  -—Du^{x^) 
ia 


J_ 

Re 


(D  ^  -  a^)  +  i(o)  -  at/. ) 


M '(X^) 


dx^ 


where 


(3.52) 


where  Dj  is  the  Chebyshev  collocation  derivative  in  computational  coordinates.  Note  that 
the  transformation  between  Cartesian  and  computational  coordinates  is  one-to-one  in  the 
wall-normal  direction  due  to  non-dimensionalization.  After  the  velocity  and  pressure 
mode  shapes  were  computed,  they  were  normalized  by  the  maximum  value  of  the  modulus 
of  . 

Multiple  disturbance  modes,  with  wavenumbers  of  a  ,  2a  ,  3a  ,  etc.,  were  obtained 
using  this  procedure  for  each  combination  of  Reynolds  number  and  a  desired.  The  initial 


72 


condition  was  computed  by  scaling  each  disturbance  mode  by  the  desired  amplitude,  e  , 
and  summing  the  base  flow,  given  by  Eq.  (2.20),  and  the  disturbances. 


^max  r 

u\xl,Xj^)  =  U^{xJ)  +  ^Realje^  ulixj)  e'*"""' 

k  =  \ 


3/^^!  Jikax! 

k=\ 
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P  I  x.  V  j 

■'  2 


u^{x.,Xj)  =  «/(-^/)  e' 


X  1 1  /  3  1  \koLx- 

=  L  Realj  ^  ,  e 

k=l  [  V  ^  2) 


(3.53) 


The  velocities  were  then  converted  from  Cartesian  to  contravariant  components. 

=  \[^^u\xl,xj) 


(3.54) 


3.4.2.  Discrete  Linear  Stability  Modes  for  Curved  Channel  Flow 

Similarly,  for  the  curved  channel,  the  mode  shapes  are  computed  using: 


a\r)  =  -—DirWir)) 
ia 
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ira 
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(3.55) 


and  D  =  —  = 

dr  ^ 
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The  mode  shapes  were  normalized  by  the  maximum  value  of  the  modulus  of  m  ®( r)  .  The 
initial  condition  was  then  computed  by  adding  the  base  flow,  given  by  Eq.  (2.27),  and  the 
disturbances.  After  the  mode  shapes  are  computed,  the  velocities  are  transformed  from 
polar  components  to  contravariant  components.  These  two  steps  are: 

k. 

'Wnax  r  V 

=  f/g(r  )  +  I^Realje^  M®(r)  e'^“®'j 

k=\  ^ 

k 

^max  /  -o'! 

=  U^(r  )  +  j  (3-56) 
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k 

max  / 

P  (6, >04)  =  Pki^j^l) 


u\d.,rj) 


(3.57) 


3.5.  Solution  Output 

The  solution  to  Eq.  (3.49)  consists  of  the  Fourier  transformed  values  of  the 
velocity  components  and  the  pressure  at  all  grid  points  for  each  time  step.  First,  the 
inverse  Fourier  transform  is  used  to  obtain  the  solution  in  physical  space.  This  solution  is 
then  output  in  several  forms  so  that  the  results  can  be  studied  effectively.  Of  course, 
planes  of  data  at  certain  times  are  output  so  that  instantaneous  velocity  vectors  and 
contours,  and  pressure  contours  may  be  shown.  For  these  results,  pressure  is  spectrally 
interpolated  from  the  half  points  to  the  grid  points.  In  these  results,  the  perturbation 
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quantities  are  also  output  by  subtracting  the  original  base  flow  from  the  instantaneous 
results.  The  evolution  of  the  disturbance  as  it  grows  or  decays  and  possibly  changes  form 
can  be  investigated  by  studying  snapshots  of  the  flow  at  various  times.  In  addition,  the 
vorticity  is  computed  for  each  snapshot  of  the  flow.  These  computations  are  summarized 
in  the  following  section. 

It  is  desirable  to  perform  further  analysis  of  the  flow  solution  every  time  step  or 
every  few  time  steps,  so  that  the  temporal  evolution  of  the  disturbances  may  be  tracked 
closely.  This  is  done  using  spectral  analysis.  The  flow  is  decomposed  into  its  Fourier 
components  in  the  periodic  directions.  The  change  in  the  amplitude  and  phase  of  each 
component  in  time  gives  information  about  the  propagation  of  waves  through  the  domain. 
The  integrated  kinetic  energy  of  each  wave  is  also  computed  and  is  used  as  the  standard 
measure  of  the  growth  or  decay  of  each  wave  and  the  balance  of  energy  between  the 
various  modes.  The  kinetic  energy  is  also  useful  in  that  it  gives  an  estimate  of  the 
truncation  error  of  the  discretization.  In  theory,  the  integrated  kinetic  energy  of  all 
neglected  modes  should  be  less  than  the  kinetic  energy  of  the  highest  retained  mode.  The 
computation  of  the  temporal  evolution  of  waves  in  the  solution  is  outlined  in  §3.5.2. 

Finally,  it  is  often  desirable  to  interpolate  the  solution  from  a  coarser  grid  to  a  finer 
grid.  This  is  needed  when  the  resolution  needs  to  be  increased  after  some  time.  It  is  also 
useful  for  displaying  results  obtained  on  a  fairly  coarse  grid.  The  low  number  of  points 
needed  in  the  periodic  directions  hides  makes  the  solution  appear  to  be  very  coarse  when 
it  is  in  fact  quite  accurate.  The  stability  results  shown  in  this  study  are  all  interpolated  on 
much  finer  grids  in  the  streamwise  direction.  Again,  spectral  interpolation  is  used.  Since 
the  method  is  very  similar  to  the  interpolation  used  to  transfer  values  from  the  grid  points 
to  the  half  points  and  vise  versa,  as  shown  in  §3.1.5,  it  will  not  be  further  elaborated  upon. 
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3.5.1.  Summary  of  Output  Flow  Quantities 

The  first  step  in  outputting  the  flow  solution  is  to  transform  the  solution  from 
Fourier  space  to  physical  space  in  the  spanwise  direction.  The  solution  is  computed  at  the 
Fourier  collocation  points  given  by  Eq.  (3.8),  with  N  =  N2,hy  multiplying  the  solution  by 
the  transform  operator,  Eq.  (3.13).  The  second  step  is  to  interpolate  the  pressure,  in  the 
wall-normal  direction,  onto  the  grid  points  using  the  interpolation  operator  Eq.  (3.32). 
These  operations  are  summarized  below: 


V"  =  Cp  U"'  ,  m=l,2,3 


P  =  K  C  'f  P 


(3.58) 


After  transformation  from  Fourier  space,  contravariant  velocity  components  are 
available  on  the  grid  points.  To  plot  velocity  vectors,  the  velocities  need  to  be  output  as 
Cartesian  components.  These  are  calculated  using: 


^  m 

u 


ax'” 


(3.59) 


However,  to  plot  contours  of  velocity  for  the  curved  channel  case,  it  is  desirable  to 
compute  the  polar  components  of  velocity,  which  can  be  computed  from  the  contravariant 
components  by  a  linear  scaling: 


(3.60) 
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The  disturbance  velocities  and  pressure  are  also  obtained  by  subtracting  the  base  flow, 
which  is  given  in  Cartesian  velocity  components  for  the  plane  channel  and  in  polar 
components  for  the  curved  channel. 

Finally,  the  vorticity  is  computed  from  the  velocities.  Cartesian  components  of  the 
vorticity  are  needed: 


d)  ^ijk 


dUj 

dx ' 


or 


- 


dUj  du. 


(3.61) 


dx‘  dxj 


,  for  /, 7, Acyclic 


The  disturbance  vorticity  is  computed  by  subtracting  the  vorticity  of  the  base  flow  from 
the  total  vorticity. 

3.5.2.  Spectral  Analysis 

The  extraction  of  information  about  the  propagation  and  growth  of  instability 
waves  in  the  solution  is  not  a  trivial  task.  The  solution  is  already  Fourier  transformed  in 
the  spanwise  direction.  The  Fourier  transform  of  the  solution  in  the  streamwise  direction 
must  be  taken: 


(3.62) 


The  next  step  is  to  integrate  some  quantity  in  the  wall-normal  direction.  One  quantity  of 
interest  is  the  kinetic  energy.  The  kinetic  energy  of  a  single  streamwise-spanwise  Fourier 
mode  is  given  by: 
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f  ’  for  k  =  Q,N^I2,  l!=0,iV,/2 
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(3.63) 


where  the  contributions  from  the  negative  modes  have  been  summed  with  the  positive 
modes. 

The  kinetic  energy  of  a  mode  exhibits  growth  or  decay  as  the  solution  proceeds. 
Linear  theory  suggests  that  this  growth  is  exponential.  Therefore,  local  exponential 
growth  is  assumed  when  computing  the  growth  rate  of  each  mode.  Assuming  that  the 
local  growth  is  given  by: 

£..('2)  (3.64) 


where  t ,  and  1 2  are  two  arbitrary  times,  and  (o  ,  is  the  imaginary  part  of  the  complex 
growth  rate,  0)^(1  may  be  computed  using: 
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(3.65) 


Obviously,  if  the  growth  is  not  actually  exponential,  the  above  equation  is  only  accurate  as 
I2  approaches  tj .  The  kinetic  energy  and  growth  rate  of  each  mode  is  computed  and 
output  at  regular  intervals. 

Because  the  kinetic  energy  is  computed  from  the  modulus  of  the  velocity  vector,  it 
contains  no  information  about  the  phase  of  each  mode.  Consequently,  the  rate  of 
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oscillation,  or  the  real  part  of  the  complex  growth  rate,  cannot  be  computed  from  the 
energy.  For  this  purpose,  the  transformed  wall-normal  contravariant  velocity  component 
was  integrate  in  the  wall-normal  direction.  This  quantity  is  labeled /i^j : 

1 

/w  =  k=-N^/2,N^/2  ,  e=iV/2,A^,/2  (3.66) 

-1 

Again  assuming  local  exponential  growth,  this  integrated  complex  quantity  should  be 
descriptive  of  the  wave  nature  of  the  flow  and  should  behave  according  to; 

(3.67) 


where  is  complex.  To  split  this  relation  into  expressions  for  the  real  and  imaginary 
parts  of  (o^j,  the  function  is  separated  into  amplitude  and  phase: 
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(3.68) 


The  real  and  imaginary  parts  of  the  complex  growth  rate  can  then  be  found  from  the  phase 
and  amplitude: 
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As  will  be  discussed  in  Section  4,  the  imaginary  part  of  the  growth  rate  computed 
from  the  two  methods  above  were  very  close,  but  not  identical.  However,  the  real  part  of 
the  growth  rate  was  found  to  match  the  theoretical  oscillation  rate  very  well  even  for  cases 
for  which  linear  theory  does  not  hold. 

3.6.  Summary 

This  chapter  contains  a  description  of  the  majority  of  the  effort  performed  in  this 
study.  The  selection  and  implementation  of  a  numerical  method  for  solving  the  governing 
equations  of  Section  2  was  one  of  the  primary  goals  of  this  study.  The  methodology  for 
studying  temporal  stability  problems  was  also  developed,  specifically,  the  computation  of 
the  initial  condition  from  linear  stability  theory.  Finally,  tools  were  developed  for  the 
analysis  of  the  results  to  determine  the  time  evolution  of  the  stability  waves  which  were 
placed  in  the  flow  at  the  initial  time.  The  following  chapter  shows  results  that  were  used 
to  validate  the  above  numerical  procedure  for  plane  and  curved  channel  flows  to  insure 
that  the  equations  are  being  solved  accurately  and  that  linear  stability  theory  results  are 
reproduced  for  small  amplitude  disturbances.  The  last  section  of  the  results  moves  beyond 
the  validation  level  and  investigates  some  early  nonlinear  flow  evolutions.  However,  the 
development  and  concurrent  validation  of  the  numerical  method  of  this  chapter  is  the  key 
result  of  this  work  and  will  allow  future  studies  to  progress  in  the  analysis  of  the  stability 
of  general  2-D  channel  flows. 


Section  4.  Plane  and  Curved  Channel  Results 

The  first  phase  of  this  study  was  the  development  and  validation  of  a  two- 
dimensional  computer  code  for  performing  the  numerical  simulation  of  the  temporal 
stability  of  channel  flows.  The  mathematical  formulation  and  numerical  approach 
described  in  the  previous  chapters  were  employed.  This  code  was  studied  extensively  to 
examine  its  effectiveness  for  simulation  of  channel  flow  problems.  Plane  and  curved 
channel  flows,  described  in  Section  2,  were  examined. 

To  generate  a  two-dimensional  version  of  the  formulation,  the  flow  solution  was 
assumed  to  be  constant  in  the  spanwise,  or  ^  ^-direction,  and  all  terms  representing 
differentiation  in  the  spanwise  direction  were  dropped.  In  the  2-D  stability  analysis,  the 
Orr-Sommerfeld  equation,  Eqs.  (2.33)  or  (2.36)  were  used  to  provide  Tollmien- 
Schlichting  wave  disturbances  for  the  plane  and  curved  channel,  respectively.  Since  the 
spanwise  terms  were  neglected,  the  temporal  stability  behavior  only  reflects  instability  in 
the  plane  of  the  streamwise  (^  *)  and  wall-normal  (^  ^)  directions.  The  important  effects  of 
spanwise  linear  stability  modes  and  the  spanwise  secondary  instability  of  Tollmien- 
Schlichting  waves  are  not  considered  in  this  study.  In  the  numerical  formulation,  requiring 
the  solution  to  be  two-dimensional  is  equivalent  to  using  a  single  Fourier  mode,  the  k=0 
mode,  to  represent  the  spanwise  variation  of  the  solution.  In  this  chapter,  the  flow 
solution  quantities  are  assumed  to  correspond  to  the  k=0  mode,  and  both  the  k  subscript 
and  ^  notation  are  dropped,  i.e.,  * . 
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First,  the  analytical  test  function  analysis  described  in  §2.1.6  was  employed  to 
determine  the  truncation  error  behavior  as  the  grid  resolution  and  time  step  were  refined. 
Spectral  accuracy  was  obtained  as  the  number  of  grid  points  in  the  streamwise  and  wall- 
normal  directions  was  increased.  Second-order  accuracy  was  obtained  as  the  time-step 
was  refined.  This  analysis,  is  shown  for  the  plane  channel  for  two  Reynolds  numbers  in 
§4.1. 

As  a  second  validation  of  the  code,  and  to  determine  the  ability  of  the  code  to 
solve  stability  problems,  linear  stability  results  were  obtained  for  both  the  plane  and  curved 
channel.  For  the  plane  channel,  two  Reynolds  numbers  were  selected  to  provide  both  a 
linearly  unstable  and  a  linearly  stable  test  case.  For  the  curved  channel,  two  cases  with 
different  radii  of  curvature  were  selected,  again  to  provide  both  a  linearly  unstable  and  a 
linearly  stable  test  case.  For  these  four  cases,  simulation  results  were  obtained  for  ten  time 
steps,  and  the  growth  and  oscillation  rates  of  the  stability  mode  were  computed.  These 
were  compared  to  linear  theory.  By  varying  the  grid  resolution,  time  step,  and  amplitude 
of  the  initial  disturbance,  the  appropriate  bounds  for  obtaining  growth  rates  consistent 
with  linear  theory  were  obtained.  These  results  are  shown  in  §4.2  for  the  plane  channel 
and  in  §4.3  for  the  curved  channel. 

The  final  set  of  results  was  obtained  to  study  the  nonlinear  stability  of  channel 
flows.  Nonlinear  stability  results  were  obtained  for  the  same  four  test  cases  simply  by 
increasing  the  amplitude  of  the  initial  disturbance.  In  addition,  higher  frequency  modes 
were  added  so  that  the  nonlinear  wave  interactions  between  the  modes  could  be  examined 
for  different  initial  amplitudes  of  the  harmonic  modes.  An  attempt  was  made  to  determine 
the  long-time  evolution  of  these  states.  Due  to  the  extremely  large  number  of  time  steps 
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required,  the  solution  for  each  case  was  only  integrated  in  time  for  up  to  thirty  periods  of 
the  Tollmien-Schlichting  wave.  These  results  are  shown  in  §4.4. 

4.1.  Analytical  Function  Testing 

The  first  validation  of  the  numerical  method  was  performed  using  the  forcing 
function  Q  introduced  in  §2.1.6.  This  method  allows  the  prescription  of  arbitrary 
analytical  functions,  U^^  and  which  are  enforced  to  be  solutions  of  the  Navier- 
Stokes  equations  by  the  proper  choice  of  the  forcing  function.  Using  this  method,  the 
spatial  discretization  error  in  the  streamwise  and  spanwise  directions,  as  well  as  the 
temporal  discretization  error,  may  be  analyzed.  Thus  without  performing  the  simulation  of 
a  flow  problem,  the  accuracy  properties  of  the  solution  method  can  be  scrutinized. 

Two  types  of  analyses  were  performed.  By  choosing  a  steady  state  test  solution 
the  spatial  discretization  was  investigated  without  running  even  a  single  time  step.  The 
test  solution  was  simply  substituted  into  the  discrete  system,  Eq.  (3.50),  and  the  residual 
was  computed.  This  residual  is  exactly  zero  in  the  absence  of  truncation  error;  when  non¬ 
zero,  the  residual  of  the  momentum  equation  is  equal  to  the  combined  error  in  the 
discretization  of  the  convective  (nonlinear)  term  and  the  diffusive  and  pressure  gradient 
(linear)  terms.  Since  the  test  velocity  was  always  chosen  to  be  a  divergence  free  vector 
function,  the  residual  of  the  continuity  equation  is  simply  the  divergence  of  the  velocity 
created  by  discretization  error. 

The  second  type  of  analysis  prescribed  a  solution  which  varied  in  time.  The  test 
solution  at  t  =  0  was  used  as  the  initial  condition.  Then,  the  simulation  was  performed  to  a 
time  oft  =  0.01  (using  At  =  0.001  for  the  spatial  accuracy  tests).  The  numerical  solution 
was  compared  to  the  test  solution  at  t  =  0.01.  This  analysis  provided  a  test  of  all  of  the 
steps  in  the  simulation  procedure.  For  the  most  part,  the  maximum  error  in  the  solution  at 
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t  =  0.01  matched  the  residual  from  the  steady-state  analysis  quite  well.  However,  the 
time-dependent  analysis  revealed  to  a  greater  extent  the  growth  of  round-off  error  in  the 
solution  procedure,  which  is  due  to  the  high  condition  number  of  the  system  for  small  At 
and  large  N.  This  issue  is  discussed  further  in  §4. 1 .5. 

4.1.1.  Test  Solutions 

To  determine  the  truncation  error  in  both  the  spatial  and  temporal  discretization,  a 
divergence-free  analytical  solution  was  selected  by  prescribing  the  stream  function  of  the 
solution.  The  general  form  for  the  stream  function  was  selected  as: 


l+ajSin(/:,^')+Z>jCos(7rsin(^'))|  =*= 

'Ua,ay^b,sin(k,e)\  * 

[l  +a^t+b^exp{kfi'\ 


(4.1) 


where  the  constants  tz,,  a^,  a,,  b^,  b^,  b,,  k^^,  k^,  and  k,  determine  specific  solutions.  The 
flow  solution  was  then  prescribed  as: 

fgde 

(4.2) 

u'(5',5’) 

This  general  solution  was  selected  so  that  it  would  contain  spatial  components 
(with  coefficients  aj  and  aj)  which  could  be  resolved  exactly  with  a  small  number  of 
modes.  This  required  the  use  of  a  sinusoidal  solution  in  the  streamwise  direction  and  a 
polynomial  solution  in  the  wall-normal  direction.  In  addition,  components  of  the  solution 
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(with  coefficients  Z?,  and  were  selected  which  could  not  be  resolved  exactly  and  would 
require  a  large  number  of  modes  to  achieve  machine  accuracy.  The  forcing  function 
vector  Q  required  to  make  Eq.  (4.2)  a  solution  to  the  momentum  equation  is  obtained  as; 

e  =  ^  tVp  -±V^«  (4.3) 

where  u  and  p  are  given  by  Eq.  (4.2).  This  substitution  was  performed  separately  for  the 
Cartesian  and  polar  momentum  equations  using  Mathematica  software  (Wolfram,  1991). 
The  resulting  expressions  were  stored  in  the  code  and  evaluated  for  each  point  every  time 
step.  Results  are  shown  only  for  the  plane  channel  case,  since  the  curved  channel  results 
were  very  similar. 

4.1.2.  Streamwise  Accuracy 

To  determine  the  truncation  error  in  the  Fourier  collocation  method  used  in  the 
streamwise  discretization,  two  results  were  calculated,  cases  A  and  B  in  Tab.  4.1.  Case  A 
was  obtained  by  specifying  the  steady-state  solution  given  by  Eq.  (4.2)  with  coefficients 
given  in  Tab.  4.1.  This  solution  is  shown  in  Fig.  4.1.  The  residual  was  calculated  for 
Reynolds  numbers  of  5,000  and  1,000  for  a  series  of  values  of  and  is  shown  in  Fig.  4.2. 
In  this  figure,  the  residual  varies  almost  exponentially  with  the  number  of  grid  points  for 
A^i  greater  than  approximately  24  (shown  by  a  near-linear  relationship  on  a  semi-log  plot), 
until  round-off  error  becomes  dominant,  which  occurs  for  Ny  around  48  for  the 
momentum  equations.  This  type  of  truncation  error  behavior  is  characteristic  of  spectral 
methods.  It  is  important  to  note  that  a  truncation  error  of  the  order  of  10'*®  can  be 
achieved  quite  feasibly. 


The  truncation  error  was  independent  of  the  Reynolds  number  for  the  cases 
shown.  This  is  surprising  since  the  Reynolds  number  significantly  changes  the  nature  of 
the  governing  equations.  However,  it  indicates  that  the  predominate  source  of  error  is  the 
nonlinear  term,  which  is  not  scaled  by  the  Reynolds  number.  It  is  possible  that  more 
extreme  values  of  Reynolds  number,  e.g.  Re  =  1  or  10  ® ,  would  show  a  difference  in  the 
residual,  but  this  was  not  investigated  here.  Also  note  that  the  residual  of  the  continuity 
equation  decreased  more  quickly  than  the  momentum  equations.  This  is  probably  due  to 
the  linearity  of  the  continuity  equation  and  is  consistent  with  the  conclusion  that  the 
nonlinear  terms  in  the  momentum  equation  are  the  primary  source  of  error  in  the  system. 


Case 

Re 

k, 

^3 

^3 

bt 

k, 

A 

1000 

5000 

1 

1 

1 

1 

0 

3 

0 

0 

1 

B 

5000 

1 

1 

1 

1 

0 

3 

1 

0 

1 

C 

5000 

1 

1 

1 

0 

1 

3 

0 

0 

1 

D 

5000 

1 

1 

1 

0 

1 

3 

1 

0 

1 

E 

5000 

1 

1 

1 

0 

0 

3 

1 

1 

1 

Table  4.1.  Test  Cases  for  Analytical  Function  Tests.  For  each  case,  the  values  of  the 
constants  in  Eq.  (4.1)  are  given  in  the  table. 


Max  Resld 
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Figure  4.1.  Contours  of  the  Prescribed  Solution  for  Streamwise  Discretization  Analysis 
of  Plane  Channel  Flow,  Case  A. 


Figure  4.2.  Demonstration  of  Spectral  Accuracy  in  the  Streamwise  Discretization  for 
Plane  Channel,  Case  A. 


L 
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The  solution  for  Case  B  increases  linearly  with  time  from  the  initial  state  and  can 
therefore  be  integrated  in  time  without  truncation  error.  For  this  case,  the  error  in  the 
solution  computed  after  10  time  steps,  with  At=0.001,  was  calculated.  These  results  are 
shown  in  Fig.  4.3.  They  are  qualitatively  the  same  as  in  Fig.  4.2.  That  is,  the  truncation 
error  of  the  solution  achieves  spectral  accuracy  as  the  number  of  grid  points  is  increased. 
It  is  useful  to  verify  that  the  residpal  of  the  numerical  discretization  is  a  good  indicator  of 
the  error  in  the  computed  solution,  since  the  maximum  value  of  the  residual  is  needed  to 
determine  convergence  at  each  time  step.  This  result  shows  the  direct  correspondence. 
When  the  residual  is  below  10''° ,  the  error  in  the  solution  is  well  below  10'^ .  As  an 
example  of  the  spatial  structure  of  the  error  in  the  solution,  the  error  is  shown  in  Fig.  4.4 
for  a  fairly  coarse  grid  of  N^=12,  A^3=12.  The  oscillatory  error  in  the  x  '-direction  shows 
that  the  error  in  the  solution  is  associated  with  missing  higher-order  terms  in  the  Fourier 
representation  of  the  solution. 


Figure  4.3.  Demonstration  of  Spectral  Accuracy  in  the  Streamwise  Discretization  for 
Plane  Channel,  Case  B.  The  maximum  error  in  the  solution  quantities  at  time 
t=0.01  are  shown. 
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Figure  4.4.  Contours  of  Error  in  Solution  at  f=0.01  for  Plane  Channel,  Case  B,  with 

Ni=l2,  N^=12.  The  error  consists  of  higher-order  Fourier  modes  neglected  in  the 
discretization. 


4.1.3.  Wall-Normal  Accuracy 

To  determine  the  truncation  error  in  the  Chebyshev  collocation  method  used  in  the 
wall-normal  discretization,  a  similar  analysis  was  performed.  The  steady-state  result  was 
obtained  by  specifying  the  Case  C  solution  of  Tab.  4.1.  Contours  of  this  solution  are 
shown  in  Fig.  4.5.  The  numerical  residual  was  calculated  for  a  series  of  values  of  Ny 
Fig.  4.6  shows  that  spectral  accuracy  behavior  is  observed  again  above  of  about  24  and 
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decreases  until  round-off  error  becomes  dominant,  around  A/3=44.  Increasing  the  number 
of  grid  points  above  this  level  begins  to  introduce  increased  round-off  error. 

Again  the  Reynolds  number  had  no  noticeable  effect  on  the  truncation  error. 


Figure  4.5.  Prescribed  Solution  for  Wall-Normal  Discretization  Analysis  of  Plane  Channel 
Flow,  Case  C. 
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Figure  4.6.  Spectral  Accuracy  is  Obtained  in  the  Wall-Normal  Discretization  for  Plane 
Channel,  Case  C. 


Again,  a  second  result  was  obtained  by  specifying  a  time-dependent  solution, 
shown  as  Case  D  in  of  Tab.  4.1.  The  error  in  the  solution  computed  after  10  time  steps, 
with  At=0.001,  is  shown  in  Fig.  4.7.  The  error  in  the  solution  is  shown  in  Fig.  4.8  and 
demonstrates  that  discretization  error  in  the  wall-normal  direction  is  composed  of  higher- 
order  polynomial  modes.  The  increased  round-off  error  in  the  residual  for  above  44 
results  in  a  slight  increase  in  the  error  of  the  solution. 


Figure  4.7.  Solution  Error  Due  to  Wall-Normal  Discretization  for  Plane  Channel,  Case  D. 


91 


Cases  C  and  D  demonstrate  that  the  Chebyshev  discretization  has  similar 
truncation  error  properties  as  the  Fourier  discretization.  One  question  that  these  cases 
helps  to  answer  is  whether  or  not  the  interpolation  in  the  wall-normal  direction  has  an 
adverse  effect  on  the  error.  The  interpolation  reduces  the  accuracy  of  the  continuity 
equation  and  the  pressure,  both  of  which  are  staggered  with  respect  to  the  momentum 
equation  and  velocity  field.  The  continuity  equation  actually  had  much  lower  residual  than 
the  momentum  equations  as  in  the  streamwise  case,  and  this  is  probably  due  to  the  fact 
that  the  continuity  equation  is  linear.  The  pressure  error  is  about  the  same  as  the  error  in 
the  streamwise  velocity.  Therefore,  it  seems  that  the  interpolation  has  no  noticeable  effect 
on  the  wall-normal  accuracy. 


Figure  4.8.  Contours  of  Error  in  Solution  for  Plane  Channel,  Case  D  with  A^i=12,  A^3=12. 
In  this  case  the  discretization  error  consists  of  higher-order  Chebyshev  polynomials 
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4.1.4.  Temporal  Accuracy 

To  assess  the  accuracy  in  the  temporal  discretization,  a  solution  with  exponential 
temporal  growth  was  used.  This  is  Case  E.  The  error  in  the  solution  was  computed  at  a 
time  of  t=0.01  for  a  family  of  time  step  values.  Naturally,  the  number  of  time  steps 
required  for  each  test  was  0.01/Ar.  The  grid  size  used  was  A, =12,  A^3=12.  The  prescribed 
solution  at  r=0.01  is  shown  in  Fig.  4.9.  The  solution  error  at  t=0.01  for  various  values  of 
At  is  shown  in  Fig.  4. 10.  The  results  for  error  in  velocity  show  that  the  temporal 
discretization  is  second-order  accurate,  as  it  should  be;  the  error  is  linear  with  a  slope  of 
approximately  2  on  a  log-log  plot  (two  decades  of  error  reduction  per  decade  of  time  step 
reduction)  down  to  a  time  step  of  10 ,  at  which  the  error  was  of  the  order  of  10'“’. 

The  analysis  yielded  erroneous  results  for  pressure  for  this  case,  due  to  an 
accumulation  in  the  pressure  error  in  the  implementation  of  the  forcing  function  analysis. 
Therefore  the  pressure  is  not  shown.  This  error  was  noticeable  only  when  a  large  number 
of  time  steps  were  used,  as  in  the  smaller  At  mns  in  Fig.  4.10.  However,  the  pressure  did 
not  show  this  behavior  in  other  tests,  so  it  is  assumed  that  an  inconsistency  in  the 
treatment  of  pressure  in  the  forcing  function  analysis  caused  the  error  accumulation.  In 
any  case,  this  error  did  not  affect  the  velocity  errors,  possibly  because  only  the  gradient  of 
the  pressure  is  used  in  the  calculations. 
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Figure  4.9.  Demonstration  of  Second-Order  Accuracy  in  Time  for  Plane  Channel,  Case  E. 


Figure  4.10.  Demonstration  of  Second-Order  Accuracy  in  Time  for  Plane  Channel, 
Case  E. 


The  temporal  error  is  distributed  spatially  over  the  entire  domain.  To  demonstrate 
this,  contours  of  the  error  in  the  solution  are  shown  for  At=0.001  at  t=0.01  in  Fig.  4.1 1. 
This  shows  that  even  if  high  accuracy  is  obtained  in  the  spatial  discretization,  the  solution 
will  contain  large  scale  spatially-distributed  temporal  errors. 

This  analysis  suggests  that  the  second-order  accurate  time  derivative  will  be 
sufficient  for  stability  simulations,  but  could  be  improved.  A  reasonable  time  step  for 
obtaining  a  solution  to  a  flow  problem  is  of  the  order  of  At  =  0.001,  and  many  of  the 


Figure  4.11.  Contours  of  Solution  Error  at  r=0.01  for  Plane  Channel,  Case  E.  This  error 
is  caused  solely  by  the  second-order  temporal  discretization,  and  demonstrates 
how  temporal  errors  are  distributed  spatially. 

simulations  performed  in  this  study  used  At  of  0.01  to  0.025.  This  study  shows  that  for 
the  temporal  error  to  be  as  low  as  the  spatial  discretization  error,  a  time  step  on  the  order 
of  At  =  10  should  be  used.  However,  this  proved  to  require  an  impractically  large 
number  of  time  steps.  With  a  higher-order  method  the  error  properties  would  probably 
improve  to  the  level  at  which  the  temporal  errors  matched  the  spatial  errors  for  reasonable 
values  of  At.  A  higher-order  method  could  be  implemented  by  replacing  the  second-order 
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trapezoidal  integration  with  a  fourth-order  implicit  integration,  but  it  would  require  the 
storage  of  a  total  of  five  levels.  This  modification  was  not  implemented  in  this  study,  but 
is  recommended  for  future  work.  Though  the  temporal  error  is  quite  small,  the  fact  that  it 
is  so  larger  than  the  spatial  error  has  some  small  but  noticeable  effects  on  the  stability 
results  later  in  this  chapter. 

4.1.5.  Condition  Number  Issues 

As  was  seen  in  the  above  studies,  the  accuracy  in  the  spatial  discretization  can  be 
reduced  to  the  level  of  machine  accuracy  (64-bit  precision)  for  a  relatively  small  number  of 
grid  points,  approximately  A^=44  for  these  test  cases,  compared  to  other  numerical 
discretization  schemes.  However,  the  spectral  collocation  derivative  matrices  used  in  this 
method  are  ill-conditioned,  and  the  condition  number  increases  with  (Canute  et  al. 
1988).  In  addition,  in  the  temporal  discretization,  all  terms  except  the  time  derivative  term 
are  multiplied  by  1/At.  Thus  the  condition  number  of  the  temporally  discrete  system  also 
increases  with  1/At.  In  the  remainder  of  this  study,  the  number  of  points  used  does  not 
exceed  N|=8,  iV3=64,  and  the  smallest  time  step  used  for  the  cases  studied  was  of  the  order 
of  10'^ .  The  previous  results  suggest  that  for  these  values,  the  ill-conditioning  of  the 
discrete  system  will  not  adversely  affect  the  solutions  obtained.  However,  in  future  phases 
of  this  study,  more  than  64  points  may  be  needed  in  the  wall-normal  direction,  as  highly 
nonlinear  interactions  lead  to  breakdown  of  the  flow  structures.  It  would  be  desirable  to 
remedy  the  round-off  error  growth  problems  to  a  level  at  which  128  points  could  be  used 
in  the  wall-normal  direction. 

The  condition  number  of  the  implicit  coefficient  matrix,  B,  in  Eq.  (3.40)  simply 
reflects  the  stiffness  of  the  spectral  discretization.  A  new  method  of  solving  the  same 
discrete  system  would  have  the  same  difficulties  with  ill-conditioning.  The  only  solution  is 
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to  find  ways  to  reduce  the  round-off  error  accumulation  in  the  solution  through  the  use  of 
higher  precision  computations.  It  is  possible  to  reduce  round-off  error  in  the  computation 
of  the  coefficient  matrix  and  in  obtaining  the  linear  algebra  solution  at  each  iteration  by 
rearranging  terms  and  operations  to  minimize  loss  of  precision.  However,  once  this  has 
been  done,  higher  precision  arithmetic  must  be  used  in  key  areas  in  which  precision  is 
being  lost.  The  key  sources  of  round-off  error  should  be  sought  out,  and  selective  use  of 
higher  precision  arithmetic  should  be  used  when  necessary  to  reduce  round-off  error 
accumulation. 

4.2.  Comparison  to  Linear  Stability  Theory  for  Plane  Channel  Flow 

The  analysis  in  §4.1  served  to  validate  and  examine  the  numerical  formulation  and 
its  implementation  in  the  two-dimensional  code.  Subsequently,  the  ability  of  the  code  to 
reproduce  linear  stability  results  for  the  plane  channel  case  was  examined.  A  Tollmien- 
Schlichting  wave  was  generated  for  the  initial  condition  as  described  in  §3.4.1.  The 
amplitude  of  the  wave  was  selected  as  A,=10  '^.  For  this  small  amplitude,  it  was  expected 
that  linear  stability  theory  would  correctly  predict  the  time  evolution  of  the  wave  as  the 
solution  was  integrated  forward  in  time  from  the  initial  condition.  Two  plane  channel 
cases  were  considered  and  are  shown  as  Case  1  and  Case  2  in  Tab.  4.2.  Case  1,  with  a 
streamwise  wavenumber  of  a=1.0  and  a  Reynolds  number  of  /?c=8,000,  has  a  small 
positive  growth  rate  and  is  therefore  slightly  unstable.  Case  2,  on  the  other  hand,  with 
a=1.25  and  /?e=4,000  has  a  negative  growth  rate  and  is  stable. 
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Case 

Re 

a 

Rc 

(Or 

(0,x  1,000 

T 

1 

8,000 

1.0 

plane 

0.247075 

2.66417 

25.4302 

2 

4,000 

1.25 

plane 

0.381718 

-10.9228 

16.4603 

3 

5,000 

100 

100 

0.374757 

3.28679 

16.7660 

4 

4,000 

5,000 

10,00 

0.530779 

-16.5727 

11.8376 

Table  4.2.  Temporal  Stability  Simulation  Test  Cases.  The  linear  stability  results  were 
obtained  from  Herbert  (1988a)  and  gives  the  complex  eigenvalue  co=(Or+/Oi),, 
where  (Or  is  the  temporal  frequency  and  (Oj  is  the  temporal  growth  rate.  The 
period  T  is  27t/o)r. 


The  solution  was  then  obtained  for  a  single  period  of  the  Tollmien-Schlichting 
wave  and  compared  with  linear  stability  theory.  In  addition,  the  grid  resolution  in  the 
wall-normal  direction,  the  time  step,  and  the  amplitude  of  the  initial  disturbance  were 
varied  parametrically,  and  the  results  were  compared  with  linear  stability  theory. 

4.2.1.  Linear  Stability  Theory  Results 

The  Tollmien-Schlichting  wave  is  an  eigenmode  of  the  linearized  governing 
equations  expanded  in  a  Fourier  series  in  the  streamwise  direction.  For  a  particular  set  of 
flow  parameters  Re  and  a,  a  finite  number  of  discrete  eigenmodes  and  eigenvalues  of  the 
linearized  system  exist.  The  eigenmode  with  the  least  stable  eigenvalue  (with  the 
maximum  imaginary  part)  is  selected  as  the  primary  mode.  Using  the  linear  stability  code 
of  Herbert  (1988a)  with  A^3=64,  the  most  stable  eigenmode  was  determined  for  Case  1  and 
Case  2.  The  resulting  eigenvalues  are  shown  in  Tab.  4.2.  The  eigenvalue  (0=(0r  -i-  /to, 
consists  of  the  temporal  frequency,  Wr  ,  and  the  growth  rate  O), .  The  complex  mode 
shapes  are  shown  in  Fig.  4.12  and  Fig.  4.13  for  Case  1  and  Case  2,  respectively.  To 
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obtain  the  Tollmien-Schlichting  wave  initial  condition,  the  mode  profiles  were  used  in 
conjunction  with  Eq.  (3.53).  Contours  of  Cartesian  velocity  components  and  pressure  are 
shown  in  Fig.  4.14  for  Case  1  and  Fig.  4.15  for  Case  2.  All  solutions  shown  in  this 
chapter  are  for  disturbance  quantities  only,  unless  otherwise  stated.  The  solution  was 
computed,  using  A^3=64,  Ar=10  '^  T,  for  a  single  period  of  the  Tollmien-Schlichting  wave. 
Results  are  shown  at  t=0.17’  and  t=T  for  Case  1  in  Figs.  4.16  and  4.17.  The  kinetic  energy 
of  the  disturbance  was  computed  using  Eq.  (3.63).  In  Fig.  4.18  this  is  compared  with  the 
growth  predicted  by  linear  theory,  which  is  given  by: 


_  2o,  ( 

~E^) 


(4.4) 


These  results  show  the  2-D  structure  of  Tollmien-Schlichting  waves  in  plane 
channel  flow.  As  time  progresses,  the  entire  wave  solution  convects  through  the  domain. 
After  one  period,  the  disturbance  contours  look  exactly  the  same  as  the  initial  Tollmien- 
Schlichting  wave  disturbance.  However,  the  amplitude  of  the  wave  has  grown  a  small 
amount  as  is  evident  in  the  increase  of  kinetic  energy.  The  solution  closely  matches  the 
predictions  of  linear  theory,  as  can  be  seen  in  Fig.  4.18.  The  slight  error  in  the  numerical 
solution  can  be  resolved  by  using  a  smaller  time  step;  however,  a  time  step  of  Ar=10'^  J  is 
already  being  used,  and  using  a  lower  time  step  would  be  prohibitively  expensive  in 
computer  time.  This  problem  is  was  also  observed  in  the  forcing  function  analysis  of  §4. 1. 
To  improve  the  accuracy  of  the  temporal  resolution  without  reducing  the  time  step,  a 
higher-order  discretization  in  time  should  be  used. 
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Figure  4.12.  Mode  Shapes  for  Plane  Channel,  Case  1. 
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Figure  4.13.  Mode  Shapes  for  Plane  Channel,  Case  2. 
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Contours  of 


Contours  of  u® 


Figure  4.14.  Contours  of  Disturbance  Quantities  for  Plane  Channel,  Case  1. 
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Figure  4.15.  Contours  of  Disturbance  Quantities  for  Plane  Channel,  Case  2. 
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Contours  of 


Contours  of  u® 


Contours  of  p 


Figure  4.16.  Contours  of  Disturbance  for  Plane  Channel,  Case  1,  t=0.\T. 
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Figure  4.17.  Contours  of  Disturbance  for  Plane  Channel,  Case  1,  r=l  .Or. 
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Figure  4.18.  Growth  of  Disturbance  Mode  Kinetic  Energy  Compared  to  Linear  theory, 
Plane  Channel,  Case  1. 
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4.2.2.  Temporal  Analysis  Results 

The  computed  solution  was  analyzed  to  determine  the  numerical  values  of  the 
temporal  frequency  Wr  and  growth  rate  to,  using  the  method  of  §3.5.  In  this  method,  the 
solution  velocity  field  is  transformed  into  Fourier  space  in  the  streamwise  direction.  Then 
the  kinetic  energy  and  the  integrated  normal  velocity  component  are  used  to  determine  the 
growth  rate  and  rate  of  oscillation  of  each  Fourier  component.  Below,  the  growth  rate 
and  rate  of  oscillation  of  the  primary  mode  are  compared  to  linear  theory.  The  analysis 
was  repeated  for  a  family  of  values  of  N^,  with  At=0.001  and  with  amplitude  of  the 
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primary  mode,  AplO  '^.  These  results  are  shown  in  Figs.  4.19  and  4.20  and  are 
summarized  in  Tab.  4.3. 

A  similar  study  was  performed  with  varying  values  of  At.  These  results  are  shown 
in  Figs.  4.21  and  4.22  and  are  summarized  in  Tab.  4.4.  They  show  that  the  temporal 
frequency  and  growth  rate  are  not  sensitive  to  the  time  step.  The  exception  is  that  the 
growth  rate  computation  based  on  kinetic  energy  loses  precision  for  At  of  10'^  or  smaller. 
This  is  due  to  the  fact  that  the  kinetic  energy  growth  rate  is  computed  from  the  square  of  a 
small  number  and  loses  digits  of  precision. 

Finally,  a  study  was  performed  to  determine  the  effect  of  disturbance  amplitude  on 
the  predicted  linear  growth.  A  family  of  initial  disturbance  amplitudes  was  used  with 
Ar=0.001  and  ^3=48.  The  results  are  shown  in  Figs.  4.23  and  4.24  and  are  summarized  in 
Tab.  4.4.  They  demonstrate  that  amplitudes  of  0.1  or  higher  are  needed  before  the  initial 
growth  of  the  wave  deviates  from  the  prescription  of  linear  theory. 

These  results  are  significant  in  that  they  demonstrate  the  limits  of  the  numerical 
scheme.  Spatial  resolution  in  the  wall-normal  direction  is  critical;  at  least  N2,=4S  points  are 
needed.  A  time  step  of  Ar=0.01  or  smaller  should  be  used.  Finally,  a  disturbance  of 
amplitude  Ai=10‘^  or  smaller  should  be  used.  Using  very  small  amplitudes  (A, =10’^)  did 
not  seem  to  cause  any  numerical  problems.  Finally,  it  is  observed  that  the  eigenvalue 
results  from  linear  theory  can  be  matched  very  closely  with  the  present  numerical  scheme, 
when  the  simulation  proceeds  within  the  ranges  demonstrated  in  this  analysis. 
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Case  1,  /?^=8,000 

Case  2,  i?e=4,000 

N, 

Ur 

xl,000 

u,(£,) 

X  1,000 

Ur 

<0,  (f>) 

X  1,000 

u,  (£,) 

X  1,000 

24 

0.247146 

2.74013 

2.66176 

0.381416 

-10.6454 

-10.9411 

32 

0.247118 

2.66204 

2.66512 

0.381698 

-10.8633 

-10.9232 

40 

0.247080 

2.66077 

2.66448 

0.381715 

-10.9185 

-10.9228 

48 

0.247076 

2.66407 

2.66445 

0.381718 

-10.9228 

-10.9228 

56 

0.247075 

2.66449 

2.66443 

0.381718 

-10.9229 

-10.9228 

exact 

0.247075 

2.66417 

2.66417 

0.381718 

-10.9228 

-10.9228 

Table  4.3.  Effect  of  Spatial  Resolution  on  Computed  Temporal  Frequency  and 

Temporal  Growth  Rate  (o,  for  Plane  Channel,  Note  that  two  different  methods  of 
computing  the  temporal  growth  rate  were  used.  The  first  method  uses  the  integral 
of  the  normal  velocity  component,  while  the  second  uses  the  kinetic  energy.  See 
§3.5  for  details. 


Case  \,Re= 

:8,000 

Case  2,  Re^ 

4,000 

At 

^R 

«,(£,) 

a)R 

(Oj  (£^j) 

X  1,000 

X  1,000 

X  1,000 

X  1,000 

10' 

0.247063 

2.66400 

2.66431 

0.381672 

-10.9182 

-10.9175 

10-' 

0.247075 

2.66412 

2.66441 

0.381718 

-10.9228 

-10.9227 

10^ 

0.247076 

2.66407 

2.66445 

0.381718 

-10.9228 

-10.9228 

10^ 

0.247076 

2.66407 

2.66747 

0.381718 

-10.9228 

-10.9231 

exact 

0.247075 

2.66417 

2.66417 

0.381718 

-10.9228 

-10.9228 

Table  4.4.  Effect  of  Temporal  Resolution  on  Computed  Temporal  Frequency  co^  and 
Temporal  Growth  Rate  co,  for  Plane  Channel. 


Case  1,  ^e=8,000 


0),(£,) 

xl,000 


Case  2,  /?e=4,000 


<0,  (£,) 

X  1,000 


A, 


10-^ 


10* 


10  = 


10  = 


exact 


Wp 


0.247078 


0.247076 


0.247076 


0.247076 


0.247075 


X  1,000 


2.64580 


2.66399 


2.66407 


2.66407 


2.66417 


2.64174 


2.66419 


2.66441 


2.66441 


2.66417 


Wp 


0.381719 


0.381718 


0.381718 


0.381718 


0.381718 


«.(/,) 

xl,000 


-10.9531 


-10.9231 


-10.9228 


-10.9228 


-10.9228 


-10.9713 


-10.9233 


-10.9228 


-10.9228 


-10.9228 


Table  4.5.  Effect  of  Disturbance  Amplitude  on  Computed  Temporal  Frequency  (Or  and 
Temporal  Growth  Rate  to,  for  Plane  Channel. 


Figure  4.19.  Effect  of  Spatial  Resolution  on  Computed  Linear  Frequency,  Wr,  for  a)  Case 
1,  and  b)  Case  2. 


CO.X1000  ©.xIOOO 
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Figure  4.20.  Effect  of  Spatial  Resolution  on  Computed  Linear  Growth  Rate,  to,,  for  a) 
Case  1,  and  b)  Case  2.  The  solid  symbols  show  values  computed  from  /,(^')  and 
the  hollow  symbols  are  computed  from 


Figure  4.21.  Effect  of  Temporal  Resolution  on  Computed  Linear  Frequency,  (Or,  for  a) 
Case  1,  and  b)  Case  2. 
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Figure  4.22.  Effect  of  Temporal  Resolution  on  Computed  Linear  Growth  Rate,  w„  for  a) 
Case  1,  and  b)  Case  2.  The  solid  symbols  show  values  computed  from  /,(^')  and 
the  hollow  symbols  are  computed  from 


b)  Case  2,  Re=4000 


Figure  4.23.  Effect  of  Disturbance  Amplitude  on  Computed  Linear  Frequency,  Wr,  for  a) 
Case  1,  and  b)  Case  2. 
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Figure  4.24.  Effect  of  Disturbance  Amplitude  on  Computed  Linear  Growth  Rate,  tOj,  for 
a)  Case  1,  and  b)  Case  2.  The  solid  symbols  show  values  computed  from 
and  the  hollow  symbols  are  computed  from 


4.3.  Comparison  to  Linear  Stability  Theory  for  Curved  Channel 
Flow 

Linear  stability  results  were  also  obtained  for  curved  channel  flows.  Cases  3  and 
4,  summarized  in  Tab.  4.2,  were  used  to  determine  the  accuracy  of  the  linear  stability 
calculations  for  the  curved  channel  geometry.  Case  3,  with  Re=5,000  and  =100  is 
unstable  due  to  the  curvature  of  the  geometry  (the  plane  channel  case  with  the  same 
parameters  is  stable).  Case  4,  with  /?e=4,000  and  R^  =1,000  is  stable.  It  is  similar  to  Case 
2  and  will  serve  as  a  useful  demonstration  of  the  effects  of  slight  curvature.  The  results 
are  similar  to  those  for  the  plane  channel. 
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4.3.1.  Linear  Stability  Theory  Results 

The  eigenvalues  for  Cases  3  and  4  are  shown  in  Tab.  4.2.  The  polar  velocity 
component  and  pressure  mode  shapes  are  shown  in  Fig.  4.25  for  Case  3  and  in  Fig.  4.26 
for  Case  4.  The  mode  profiles  were  used  to  calculate  the  curved  channel  Tollmien- 
Schlichting  waves  using  Eq.  (3.56).  Contours  of  the  disturbance  polar  velocity 
components,  and  u' ,  along  with  the  pressure  p,  are  shown  in  Fig.  4.27  for  Case  3  and 
Fig.  4.28  for  Case  4.  Again,  the  solution  was  computed,  using  Nj=64,  At=lO'^T,  for  a 
single  period  of  the  Tollmien-Schlichting  wave.  The  kinetic  energy  of  the  disturbance 
mode  was  also  calculated  and  is  compared  to  the  kinetic  energy  predicted  by  linear  theory. 
This  is  shown  in  Fig.  4.29. 
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Figure  4.26.  Mode  Shapes  for  Curved  Channel,  Case  4. 
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Figure  4.27.  Disturbance  Quantities  for  the  Tollmien-Schlichting  Wave  Initial  Condition 
for  the  Curved  Channel,  Case  3. 
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Contours  of  u® 


Contours  of  u"^ 


Contours  of  p 


Figure  4.28.  Disturbance  Quantities  for  the  Tollmien-Schlichting  Wave  Initial  Condition 
for  the  Curved  Channel,  Case  4. 
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Figure  4.29.  Disturbance  kinetic  energy  for  one  period  T  for  Curved  Channel,  Case  3. 
Numerical  simulation  is  shown  as  dashed  hne,  linear  stability  solution  shown  as 
dotted  line. 


4.3.2.  Temporal  Analysis  Results 

The  computed  solution  was  analyzed  to  determine  the  numerical  values  of  the 
temporal  frequency  Wp.  and  growth  rate  tOj  for  the  curved  channel  as  was  done  for  the 
plane  channel  in  §4.2.2.  The  spectral  analysis  of  §3.5  was  used  to  obtain  results  for  a 
family  of  values  of  iVj  with  At=0.001,  Aj=10'^  These  results  are  shown  in  Figs.  4.30  and 
4.31  and  are  summarized  in  Tab.  4.6. 
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The  results  were  also  obtained  for  varying  values  of  At  with  A^3=48,  Ai=10‘^. 

These  results  are  shown  in  Figs.  4.32  and  4.33  and  are  summarized  in  Tab.  4.7.  Again  the 
growth  rate  to,  is  more  sensitive  to  resolution  than  the  temporal  frequency  oOr. 

Finally,  a  family  of  initial  disturbance  amplitudes  was  used  with  ^3=48,  Ar=10'^. 
The  results  are  shown  in  Figs.  4.34  and  4.35  and  are  summarized  in  Tab.  4.8. 

These  studies  confirm  that  the  stability  solutions  for  curved  channel  flow  behave 
much  the  same  as  those  for  plane  channel  flow.  However,  in  some  cases,  as  few  as  three 
four  digits  of  accuracy  are  obtained  in  either  the  temporal  frequency  and  growth  rate  when 
the  normal-velocity  integral  was  used  to  compute  the  evolution  of  the  wave.  The  results 
based  on  kinetic  energy  are  correct  to  at  least  five  decimal  places  in  comparison  to  the 
linear  theory  for  the  cases  shown.  Again  the  results  suggest  that  the  discretization 
requires  at  least  a  wall-normal  resolution  of  ^3=48  points,  a  temporal  resolution  of 
At=10'^  or  smaller,  and  an  amplitude  of  A, =10'^  or  smaller.  Note  that  in  the  curved 
channel  case,  however,  disturbances  amplitudes  of  A,=10'^  or  smaller  began  to  cause 
numerical  difficulties  for  Case  4.  It  was  observed  that  for  radii  of  curvature  even  larger 
than  =  1,000,  numerical  errors  began  to  arise  in  all  calculations.  The  metric  terms 
defining  the  curved  channel  geometry  lose  precision  for  high  R^.  To  simulate  very  large 
radii  of  curvature,  the  curved  channel  geometry  needs  to  be  reformulated  so  that  distances 
are  non-dimensionalized  with  the  channel  height,  h,  instead  of  the  radius  of  curvature,  R^ . 
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Case  3,  Re= 

5,000,  R=m 

Case  4,  /?e=4,000,  Rc=l  ,000 

N, 

(Or 

G)j  (£'j) 

(Or 

(Oi  (E^) 

X  1,000 

xl,000 

xl.OOO 

X  1,000 

24 

0.374860 

3.44599 

3.28135 

0.530855 

-16.1705 

-16.5948 

32 

0.374845 

3.29589 

3.28771 

0.530914 

-16.5031 

-16.5704 

40 

0.374790 

3.28305 

3.28685 

0.530800 

-16.5703 

-16.5726 

48 

0.374782 

3.28548 

3.28680 

0.530781 

-16.5723 

-16.5728 

56 

0.374782 

3.28550 

3.28679 

0.530779 

-16.5727 

-16.5728 

exact 

0.374757 

3.28679 

3.28679 

0.530779 

-16.5727 

-16.5727 

Table  4.6.  Effect  of  Spatial  Resolution  on  Computed  Temporal  Frequency  0)^  and 
Temporal  Growth  Rate  cO]  for  Curved  Channel. 


Case  3,  /?e=5,000,  /?,=100 

Case  4,  /?e=4,000,  Rc=  1,000 

At 

(0]  (£*1) 

^r 

(£^i) 

X  1,000 

X  1,000 

X  1,000 

X  1,000 

10’ 

0.374738 

3.28598 

3.28674 

0.530655 

-16.5602 

-16.5581 

10-' 

0.374817 

3.28564 

3.28680 

0.530779 

-16.5723 

-16.5726 

10-' 

0.374822 

3.28548 

3.28680 

0.530781 

-16.5723 

-16.5728 

10-'’ 

0.374823 

3.28546 

3.28680 

0.530781 

-16.5723 

-16.5728 

10-^ 

0.374823 

3,28546 

3.28681 

0.530781 

-16.5723 

-16.5728 

exact 

0.374757 

3.28679 

3.28679 

0.530779 

-16.5727 

-16.5727 

Table  4.7.  Effect  of  Temporal  Resolution  on  Computed  Temporal  Frequency  0)r  and 
Temporal  Growth  Rate  (Oj  for  Curved  Channel. 


Case  3,  ^e=5,000,  J?£=100 


Case  4,  ^?e=4,000,  ^^=1,000 


>4, 

(Or 

X  1,000 

w,  (£,) 

X  1,000 

(Or 

X  1,000 

«.  (Ed 

X  1,000 

10“ 

0.374785 

3.26480 

3.25982 

0.530782 

-16.6130 

-16.6415 

10' 

0.374782 

3.28537 

3.28653 

0.530781 

-16.5727 

-16.5735 

10-' 

0.374782 

3.28549 

3.28679 

0.530781 

-16.5723 

-16.5728 

10-^ 

0.374782 

3.28548 

3.28680 

0.530781 

-16.5723 

-16.5728 

10-= 

0.374782 

3.28548 

3.28681 

0.530781 

-16.5721 

-16.5727 

exact 

0.374757 

3.28679 

3.28679 

0.530779 

-16.5727 

-16.5727 

Table  4.8.  Effect  of  Initial  Amplitude  on  Computed  Temporal  Frequency  co^  and 
Temporal  Growth  Rate  «,  for  Curved  Channel. 


Figure  4.30.  Effect  of  Spatial  Discretization  on  Computed  Linear  Frequency  for  Curved 
Channel,  a)  Case  3,  and  b)  Case  4. 
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Figure  4.31.Effect  of  Spatial  Discretization  on  Computed  Linear  Growth  Rate,  tOj,  for  a) 
Case  3,  and  b)  Case  4.  The  solid  symbols  show  values  computed  from  and 
the  hollow  symbols  are  computed  from 


Figure  4.32.  Effect  of  Temporal  Discretization  on  Computed  Linear  Frequency  for 
Curved  Channel,  a)  Case  3,  and  b)  Case  4. 
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Figure  4.33.  Effect  of  Temporal  Discretization  on  Computed  Linear  Growth  Rate,  co,,  for 
a)  Case  3,  and  b)  Case  4.  The  solid  symbols  show  values  computed  from 
and  the  hollow  symbols  are  computed  from 


Figure  4.34.  Effect  of  Disturbance  Amplitude  on  Computed  Linear  Frequency  for  Curved 
Channel,  a)  Case  3,  and  b)  Case  4. 
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Figure  4.35.  Effect  of  Disturbance  Amplitude  on  Computed  Linear  Growth  Rate,  cOj,  for 
a)  Case  3,  and  b)  Case  4.  The  solid  symbols  show  values  computed  from  fi(V) 
and  the  hollow  symbols  are  computed  from 


4.4.  Nonlinear  Stability  Results  for  Plane  Channel  Flow 

In  addition  to  the  forcing  function  tests  and  linear  stability  validations,  a  limited 
number  of  runs  were  performed  to  validate  the  code  for  nonlinear  stability  calculations  and 
to  study  the  effects  of  nonlinearity  on  the  solution.  Plane  channel  results  were  obtained 
for  Cases  1  and  2  in  Tab.  4.2.  In  addition,  a  few  curved  channel  results  were  obtained,  but 
they  showed  few  qualitative  differences  from  the  plane  channel  results  and  will  not  be 
shown  here. 

The  results  show  the  effects  of  nonlinearity  on  the  early  flow  evolution  of  the 
primary  mode.  The  threshold  initial  amplitude  above  which  nonlinear  behavior  occurs  is 
explored.  In  addition,  for  Case  2,  it  appears  that  an  equilibrium  state  is  being  approached 
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in  the  long  time  evolution.  Plots  of  the  flow  field  after  twenty  periods  of  the  Tollmien- 
Schlichting  wave  show  that  the  flow  resembles  the  finite-amplitude  states  observed  in 
experiments  and  simulations. 

Two  validations  were  available  for  the  nonlinear  results  for  the  plane  channel.  The 
direct  simulation  code  of  Zang  (1984)  was  used  to  obtain  results  for  the  conditions  of 
Case  1.  These  results  compared  well  with  the  results  obtained  in  this  study.  In  addition,  a 
finite-amplitude  state  was  computed  in  Orszag  and  Patera  (1983)  for  the  conditions  of 
Case  2.  These  results  also  compare  well. 

4.4.1.  Early  Time  Results 

Results  were  obtained  for  Case  1  and  Case  2  for  varying  values  of  initial 
amplitude,  A  ^ .  The  kinetic  energy  of  the  primary  mode  is  shown  as  a  function  of  time  for 
these  cases  in  Figs.  4.36  and  4.37,  respectively.  The  results  show  several  interesting 
trends.  In  all  cases,  the  initial  evolution  closely  matches  linear  stability  theory,  which 
predicts  exponential  growth  (which  appears  as  linear  on  a  log  plot).  However,  during  the 
first  period  of  the  Tollmien-Schlichting  wave,  the  high  amplitude  cases  begin  to  deviate 
from  linear  evolution.  For  Case  1,  the  high  amplitude  result,  with  A  i  =  0.1,  is  amplified 
super-exponentially,  through  the  nonlinear  interactions  between  the  fundamental  mode  and 
its  harmonics.  In  contrast,  the  energy  for  A  i  =  0.01  grows  exponentially,  though  the 
growth  rate  is  slightly  higher  than  the  linear  stability  result. 

For  Case  2,  the  flow  is  linearly  stable,  and  the  disturbances  initially  decay 
according  to  linear  theory.  However,  about  half-way  through  the  first  period,  the  kinetic 
energy  of  the  primary  mode  for  initial  amplitudes  higher  than  A ,  =  0.01  begins  to  slow  its 
decay.  For  cases  with  initial  amplitudes  of  A  i  =  0.05  and  higher,  the  wave  begins  to  grow. 


123 


In  both  of  these  cases,  the  effect  of  nonlinear  interactions  is  to  destabilize  the 
primary  mode,  causing  it  to  increase  super-exponentially  for  the  linearly  unstable  case,  and 
to  reverse  its  growth  in  the  linearly  stable  case.  Referring  to  the  weakly  nonlinear  theory 
of  §1.1.3,  these  effects  confirm  that  the  Landau  coefficient,  {,  is  positive  for  the  plane 
channel.  In  this  context.  Case  2  is  a  case  of  subcritical  instability.  Though  the  flow  is 
linearly  stable,  it  is  unstable  to  disturbances  of  a  finite  size.  However,  note  that  the 
evolution  predicted  by  the  Landau  equation,  Eq.  (1.1),  differs  from  these  results  in  a 
significant  way.  The  Landau  equation  does  not  predict  that  the  mode  will  first  grow  or 
decay  according  to  linear  theory,  and  then  deviate  from  this  after  some  time.  Rather  it 
predicts  monotonic  growth  for  an  intial  amplitude  above  a  threshold  amplitude,  A  ^ .  In  the 
present  results,  A  g  is  approximately  0.06.  For  initial  amplitudes  above  A  ^ ,  the  disturbance 
first  decays,  and  then  grows.  This  suggests  that  other  mechanisms  are  at  work  during  the 
early  transient  which  do  not  allow  the  primary  mode  to  follow  the  evolution  equation 
obtained  from  weakly  nonlinear  stability  theory.  One  possible  mechanism  that  could 
explain  this  behavior  is  that  for  large  amplitudes,  the  least  stable  mode  selected  from  linear 
stability  theory,  may  not  be  the  only  mode  which  is  being  triggered  by  the  nonlinear 
interactions.  The  effect  of  this  is  the  modification  of  the  shape  of  the  primary  disturbance 
during  the  early  transient.  After  the  transient,  the  mode  shape  remains  fixed  in  shape  and 
then  evolves  nonlinearly. 
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Figure  4.36.  Growth  of  the  primary  disturbance  mode  for  one  period,  T,  of  the  Tollmien- 
Schlichting  wave,  for  the  plane  channel.  Case  1,  with  Re  =  8,000  and  a  =  1.0.  The 
results  for  A ,  =0.1  show  significant  nonlinearity. 
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Figure  4.37.  Growth  of  the  primary  disturbance  mode  for  one  period,  T,  of  the  Tollmien- 
Schlichting  wave,  for  the  plane  channel.  Case  2,  with  Re  =  4,000  and  a  =  1.25. 
The  results  for  A ,  >  0.01  show  significant  nonlinearity.  In  particular,  the  growth 
rate  becomes  positive  for  A  j  >  0.04. 


4.4.2.  Long  Time  Growth  of  the  Primary  Mode 

The  early  time  results  suggest  nonlinear  evolution  according  to  weakly  nonlinear 
stability  theory  after  a  transient  of  about  half  of  a  period  of  the  Tollmien-Schlichting  wave, 
The  next  question  to  investigate  is  whether  or  not  the  early  time  evolution  is  indicative  of 
the  long  time  behavior.  Fig.  4.38  shows  the  nonlinear  growth  of  the  high  amplitude  result 
from  Case  1.  The  super-exponential  growth  during  the  first  period  does  not  continue. 

The  growth  rate  decreases  as  the  disturbance  continues  to  grow.  This  is  consistent  with 
the  notion  that  weakly -nonlinear  ideas  only  apply  for  small  amplitude  disturbances.  Once 
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disturbances  grow  beyond  this  range,  higher-order  terms  would  be  needed  to  predict  the 
evolution  of  the  disturbance.  In  the  extreme,  as  the  disturbance  becomes  of  order  1,  the 
asymptotic  series  approach  completely  breaks  down,  and  an  infinite  number  of  terms 
would  be  needed  to  predict  the  fully  nonlinear  behavior. 

Fig.  4.39  shows  the  evolution  for  Case  2  for  the  first  five  periods.  The  case  with 
initial  amplitude  A  j  =  0.05  first  decays,  then  grows,  then  decays.  Eventually  it  will  decay 
to  zero  amplitude.  The  energy  for  initial  amplitude  A  i  =  0.08  is  growing  slightly  after  2 
periods.  After  some  additional  investigation,  not  shown  in  these  figures,  the  threshold 
amplitude,  A  ^  for  Case  2  was  determined  to  be  approximately  0.06.  The  flow  is  stable  to 
initial  disturbances  below  this  amplitude,  and  is  unstable  to  disturbances  above  this 
amplitude. 
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Plane  Channel,  Case  1,  Re=8,000,  a=1.0 


Figure  4.38.  Growth  of  the  primary  disturbance  mode  for  10  periods  of  the  Tollmien- 

Schlichting  wave,  for  the  plane  channel.  Case  1,  with  Re  =  8,000  and  a  =  1.0.  The 
growth  of  the  A ,  =  0.1  result  slows  significantly  after  the  first  period.  Note  that 
the  A ,  =  10  '^  and  A  ^  =  10  solutions  were  obtained  for  only  one  period. 
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Figure  4.39.  Growth  of  the  primary  disturbance  mode  for  5  periods  of  the  Tollmien- 

Schlichting  wave,  for  the  plane  channel,  Case  2,  with  Re  =  4,000  and  a  =  1 .0.  The 
critical  value  of  A ,  for  long  term  growth  is  approximately  0.06.  Initial 
disturbances  above  this  amplitude  will  not  decay. 


4.4.3.  Evolution  of  Harmonic  Modes 

Results  are  shown  for  Cases  1  and  2  to  demonstrate  the  nonlinear  behavior  of  the 
first  three  harmonic  modes.  The  modes  are  designated  in  the  manner  commonly  used  in 
the  literature.  The  mode  number  (k,^)  denotes  the  mode  with  a  streamwise  wavenumber, 
k,  and  a  span  wise  wavenumber  L  In  this  case,  the  primary  mode  is  denoted  (1,0),  and  its 
harmonics  are  (2,0),  (3,0),  (4,0).  For  Case  1,  with  A  j  =  0.001,  the  results  are  shown  in 
Fig.  4.40,  and  are  compared  with  the  results  of  Zang  (1984).  The  initial  conditions  were 
slightly  different  between  the  present  results  and  the  results  of  Zang,  and  this  leads  to  an 
interesting  observation.  The  initial  conditions  for  present  results  were  comprised  of  the 


129 


base  flow,  the  primary  disturbance,  and  harmonics  of  the  primary  disturbance  at  an 
amplitude  smaller  than  the  primary  mode.  For  the  results  of  Fig.  4.40,  the  initial  condition 
was  obtained  with  A  j  =  0.00 1 ,  A  2  =  0.000 1 ,  A  3  =  A  4=  0.  The  results  of  Zang  were 
obtained  with  A  j  =  0.001 ,  A2  =  A3  =  A4=0.  After  a  short-time  transient,  the  two  results 
match  fairly  well,  indicating  that  the  initial  amplitude  of  the  harmonics  is  irrelevant  for 
determining  the  evolution  of  the  flow.  Only  the  amplitude  of  the  primary  mode  is 
important,  as  it  will  quickly  force  the  harmonic  modes  at  certain  amplitudes  independent 
of  their  initial  amplitudes.  This  observation  was  reinforced  by  many  other  results  not 
shown  here. 

The  harmonics  are  shown  for  Case  2,  with  A ,  =  0.001,  for  twenty  periods  of  the 
Tollmien-Schlichting  wave.  Fig.  4.41  shows  that  the  primary  mode  decays  exponentially. 
The  second  mode  also  decays  exponentially.  However,  the  third  mode  stops  its  decay 
after  several  periods,  and  then  stabilizes  at  a  fixed  amplitude.  This  behavior  leaves  the 
issue  of  the  long-time  behavior  of  this  flow  a  mystery.  When  the  primary  mode  decays  to 
the  level  of  the  (3,0)  mode,  it  is  possible  that  it  will  also  stabilize,  through  interactions  with 
this  harmonic.  It  is  also  possible  that  the  nonlinear  interactions  will  cause  the  primary  to 
begin  growing  again.  This  suggests  that  in  the  case  of  subcritical  instability,  the  threshold 
amplitude  for  long-time  stability  may  be  quite  small.  Though  a  much  larger  disturbance 
was  required  in  order  to  observe  nonlinear  growth  (and  therefore  subcritical  instability)  in 
the  early  time  evolution,  it  is  possible  that  the  actual  threshold  amplitude  for  stability  is 
much  smaller. 

The  long  time  behavior  of  Case  2,  with  A  ,  =  0.1,  is  shown  in  Fig.  4.42.  The 
harmonics  in  this  case  all  behave  nonlinearly,  as  does  the  primary  mode.  After  20  periods 
of  the  Tollmien-Schlichting  wave,  the  growth  rate  decreased  significantly  from  the  initial 
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nonlinear  growth.  The  harmonics  behave  in  a  similar  manner  as  the  primary  mode.  The 
evolution  of  these  modes  seems  to  suggest  that  the  primary  mode  and  its  harmonics  are 
approaching  an  equilibrium  state.  This  state,  called  a  finite-amplitude  state,  is  one  in 
which  the  nonlinear  interactions  between  the  primary  mode  and  its  harmonics  reaches  an 
equilibrium.  This  represents  a  bifurcation  of  the  flow  to  a  new  state,  which  is  stable  to 
streamwise  disturbances  (note  that  the  secondary  instability  of  this  state  to  spanwise 
disturbances  will  cause  it  to  break  down  as  discussed  in  Section  1 .)  The  long  time 
evolution  of  subcritical  instability,  therefore,  is  to  a  stable  equilibrium,  given  a  large 
enough  disturbance.  As  was  mentioned  above,  it  is  possible  that  smaller  initial 
disturbances  may  also  eventually  lead  to  growth,  in  which  case,  the  same  equilibrium 
should  be  achieved  as  for  the  high  amplitude  case. 
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Figure  4.40.  Growth  of  the  primary  disturbance  mode  (1,0)  and  its  harmonics  (2,0),  (3,0), 
and  (4,0)  for  5  periods  of  the  Tollmien-Schlichting  wave,  for  the  plane  channel. 
Case  1.  The  symbols  are  the  results  of  Zang  (1988a)  for  similar  conditions,  but 
with  a  different  initial  amplitude  for  the  harmonic  (2,0). 
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Figure  4.41.  Growth  of  the  primary  disturbance  mode  (1,0)  and  its  harmonics  (2,0),  (3,0), 
and  (4,0)  for  20  periods  of  the  Tollmien-Schlichting  wave,  for  the  plane  channel. 
Case  2.  The  initial  amplitude  is  0.001.  The  (3,0)  modes  stops  decaying  after 
about  2  periods,  suggesting  that  after  a  long  time  nonlinear  interactions  may  force 
subsequent  growth  of  the  primary  mode. 
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Case  2,  Plane  Channel,  Re=4,000,  0=1.25,  A,=10'' 


Figure  4.42.  Growth  of  the  primary  disturbance  mode  (1,0)  and  its  harmonics  (2,0),  (3,0), 
and  (4,0)  for  20  periods  of  the  Tollmien-Schlichting  wave,  for  the  plane  channel. 
Case  2.  The  initial  amplitude  is  0. 1 .  Though  this  case  is  linearly  stable,  the 
disturbances  are  growing  to  a  large  amplitude.  This  is  an  example  of  subcritical 
instability. 


4.4.4.  Finite  Amplitude  State  for  Plane  Channel  Flow 

The  flow  field  for  Case  2,  with  an  initial  amplitude  A  j  =  0.1,  was  investigated  at 
the  time  t=20T.  The  results  were  in  good  agreement  with  the  computation  of  a  finite- 
amplitude  state  performed  by  Orszag  and  Patera  (1983).  Figs.  4.43  and  4.44,  respectively, 
show  contours  of  the  disturbance  qualities,  including  vorticity,  and  the  total  flow 
quantities.  In  addition,  Fig.  4.45  shows  vectors  and  streamlines  of  the  disturbance  velocity 
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field.  The  finite-amplitude  Tollmien-Schlichting  wave  is  somewhat  distorted  from  its 
initial  shape,  but  retains  its  basic  structure.  Compare  with  the  initial  shape  given  in 
Fig.  4.15.  The  disturbance  has  approximately  tripled  in  amplitude  from  its  initial 
amplitude,  and  thus  the  disturbance  is  sufficiently  large  to  significantly  effect  the  total 
streamwise  velocity  and  vorticity,  as  shown  in  Fig.  4.44.  A  pronounced  waviness  can  be 
seen  in  the  flow.  This  is  consistent  with  experimental  observations  of  instability  waves  in 
shear  flows.  They  are  typically  sufficiently  large  to  be  observed  in  flow  visualizations. 
The  streamlines  of  Fig.  4.45  appear  similar  to  the  finite  amplitude  state  computed  by 
Orszag  and  Patera  (1983)  for  the  same  Reynolds  number  and  streamwise  wavelength. 


Figure  4.43.  Contours  of  disturbance  quantities  after  20  periods  of  the  Tollmien- 

Schlichting  wave  for  the  plane  channel,  Case  2.  The  flow  is  approaching  a  finite- 
amplitude  state.  This  state  is  notably  different  than  the  original  Tollmien- 
Schlichting  wave  disturbance. 


Figure  4.44.  Contours  of  total  streamwise  velocity  and  vorticity  after  20  periods  of  the 
Tollmien-Schlichting  wave  for  the  plane  channel.  Case  2.  The  growth  of  the 
disturbance  has  significantly  altered  the  appearance  of  the  flow. 
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Figure  4.45.  Vectors  and  streamlines  of  the  disturbance  velocity  after  20  periods  of  the 
Tollmien-Schlichting  wave  for  the  plane  channel,  Case  2.  The  streamlines  match 
the  results  of  Orszag  and  Patera  (1983). 

4.5.  Summary 

The  results  of  this  chapter  serve  first  to  validate  the  numerical  approach.  They 
demonstrate  spectral  accuracy  both  in  the  overall  truncation  error,  as  shown  in  the 
analytical  function  analysis,  and  in  the  linear  stability  results.  However,  the  limitations  of 
the  current  scheme  are  also  revealed.  Round-off  error  is  a  potential  problem  for  a  large 
number  of  grid  points,  for  small  time  steps,  and  for  small  disturbance  amplitudes.  The 
only  solution  for  these  problems  is  to  use  a  higher  precision  level  in  computing  certain 
terms,  as  discussed  in  §4.1.5.  These  problems  have  not  affected  the  solutions  of  interest, 
however,  and  will  be  left  for  future  study. 

Secondly,  the  results  obtained  thus  far  have  investigated  some  of  the  effects  of 
nonlinearity  of  Tollmien-Schlichting  waves  in  2-D  plane  and  channel  flows.  These  results 
show  that  disturbance  amplitudes  of  the  order  of  1%  of  the  streamwise  mean  flow  or 
higher  are  needed  in  order  to  observe  nonlinear  behavior.  This  is  consistent  with  the 
general  consensus  in  the  literature  that  “small  amplitude”  disturbances  are  those  with 
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amplitude  less  than  1%,  implying  that  nonlinear  behavior  begins  to  occur  near  1%.  Also, 
the  threshold  amplitude  for  subcritical  instability  of  plane  channel  flow  at  a  Reynolds 
number  of  4,000  was  determined  to  be  approximately  6%,  though  the  actual  value  may  be 
much  lower  if  nonlinear  growth  appears  only  after  a  long  time.  For  amplitudes  above  this 
threshold,  the  disturbance  will  eventually  cause  the  flow  to  bifurcate  to  a  finite  amplitude 
state,  at  which  the  amplitude  of  the  primary  mode  is  approximately  30%  of  the  mean  flow. 
These  results  are  just  a  tiny  sample  of  the  nonlinear  results  that  can  be  obtained  using  the 
numerical  approach  developed  in  this  study.  Much  work  remains  to  be  done  to  employ 
the  numerical  tool  which  has  been  developed  to  further  investigate  the  nonlinear  evolution 
of  instabilities. 


Section  5.  Conclusions  and  Recommendations 

The  purpose  of  this  work  was  achieved  not  only  in  the  final  outcome,  represented 
by  the  plots  and  tables  of  results,  but  by  the  development  of  a  methodology  for  studying 
the  stability  of  shear  flows.  Further  development  and  validation  is  needed  to  complete  the 
3-D  version  of  the  computer  code.  In  addition,  more-complex  channel  flow  geometries 
need  to  implemented.  Both  of  these  additions  to  the  current  code  must  be  carefully 
validated.  These  are  the  next  steps  of  this  study.  Looking  beyond  these  steps,  however, 
recommendations  can  be  made  for  the  future  use  of  the  methodology  described  in  this 
study.  This  methodology  should  be  used  for  the  investigation  of  current  nonlinear  stability 
problems.  Three  major  areas  of  investigation  are  discussed  below,  along  with  the 
recommendations  for  study  using  the  3-D  generalized  channel  code. 

5.1.  Nonlinear  Wave  Interactions 

The  recent  advances  in  the  field  of  instability  and  transition  to  turbulence  have 
largely  focused  on  the  investigation  of  nonlinear  interactions  between  various  instability 
waves.  These  interactions  have  been  determined  in  experiments  to  be  critical  to  the 
understanding  of  many  prevalent  flow  configurations  in  which  different  instability 
mechanisms  compete  for  dominance.  Numerous  theoretical  studies,  usually  combined 
with  a  significant  amount  of  computation,  have  investigated  various  types  of  wave 
interactions  and  have  attempted  to  explain  experimental  observations.  Some  examples  of 
problems  in  which  wave  interactions  are  important  are  the  flow  over  a  swept  wing  (Reed 
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and  Saric,  1989),  the  secondary  instability  of  Tollmien-Schlichting  waves  in  a  boundary 
layer  (Orszag  and  Patera,  1983),  and  the  Tollmien-Schlichting  wave-Dean  vortex 
interaction  in  curved  channel  flows  {Singer,  Erlebacher  and  Zang,  1992)  Wave 
interactions  are  responsible  for  many  of  the  remaining  unanswered  questions  in  the  study 
of  transition  to  turbulence. 

Many  of  the  current  studies  of  wave  interactions  have  used  either  plane  or  curved 
channel  flows  as  the  geometry  of  interest.  This  is  due  to  the  fact  that  plane  channels 
exhibit  many  of  the  characteristics  of  plane  boundary  layers  and  can  be  used  to  investigate 
secondary  instability  and  nonlinear  breakdown  at  a  lower  cost  than  the  comparable 
boundary  layer  simulations  (Kleiser  and  Zang,  1989).  Also,  curved  channel  flow  is  one  of 
the  simplest  settings  for  the  study  of  curvature-induced  instabilities.  The  3-D  general 
channel  code  can  be  used  to  further  investigate  nonlinear  wave  interactions  in  plane  and 
curved  channels.  Results  of  other  authors  for  secondary  instability  and  Tollmien- 
Schlichting  wave-streamwise  vortex  interactions  in  plane  and  curved  channels  should  be 
reproduced  and  advanced.  In  addition,  the  generalized  coordinate  system  could  be  used  to 
provide  small  perturbations  on  the  channel  walls  which  could  alter  the  wave  interactions 
and  lead  to  a  different  route  to  turbulence.  This  area  has  been  investigated  very  little  in 
existing  literature. 

5.2.  Linear  and  Nonlinear  Stability  of  Complex  Flows 

In  flows  of  engineering  interest,  many  shear  layers  exist  with  different  spatial  scales 
and  different  stability  characteristics.  Recent  studies  have  investigated  sample  complex 
flows  in  which,  for  example,  regions  of  separated  flow  occur  in  addition  to  the  boundary 
layer  or  other  shear  layers.  Some  of  these  studies  are  flow  over  a  cylinder  (Tomboulides, 
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Triantafyllou  and  Karniadakis,  1992),  over  a  backward-facing  step  (Kaiktsis,  Kamiadakis 
and  Orszag,  1991),  and  in  the  “eddy-promoter”  channel,  that  is  a  channel  with  a  small 
cylinder  imbedded  in  it  to  cause  small  scale  separation  (Karniadakis,  Mikic  and  Patera, 
1988,  and  Karniadakis  and  Triantafyllou,  1992).  For  these  flows,  a  simple  analytical  base 
flow  is  not  available,  and  consequently,  the  linear  stability  characteristics  of  that  base  flow 
are  not  known.  Numerical  simulation  must  be  used  to  determine  the  base  flow,  and  then 
to  investigate  the  linear  stability.  Nonlinear  stability  results  then  follow.  Recent  results 
have  suggested  that  the  Tollmien-Schlichting  wave  phenomenon  is  very  dominant  even  in 
flows  which  are  significantly  different  from  canonical  flows.  However,  the  growth  rate  of 
the  Tollmien-Schlichting  wave  is  significantly  influenced  by  the  presence  of  additional 
shear  layers  in  the  base  flow. 

The  general  channel  geometry  is  well-suited  to  the  investigation  of  complex  flows. 
Though  the  current  spectral  method  requires  a  smooth,  analytical  coordinate 
transformation,  smooth  cavities,  grooves,  and  obstacles  can  be  placed  on  the  surface  of 
the  channel  and  can  cause  flow  separation  and  the  occurrence  of  multiple  shear  layers  of 
different  scales  within  the  flow.  A  wealth  of  information  can  be  obtained,  even  with  only  a 
2-D  simulation.  The  effects  of  spanwise  disturbances  can  only  be  investigated  after  the 
evolution  of  the  2-D  Tollmien-Schlichting  wave  is  determined.  The  first  step  in  this 
direction  would  be  to  reproduce  the  results  of  Guzman  and  Amon  (1982)  for  the  flow  in  a 
periodic  converging-diverging  channel.  A  consistent  methodology  for  computing  the  base 
flow  and  the  linear  stability  characteristics  should  be  determined.  Then,  other  geometries 
should  be  investigated.  Eventually,  3-D  simulations  should  be  performed  to  study  the 
breakdown  of  Tollmien-Schlichting  waves  due  to  spanwise  disturbances  in  these  flows. 
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5.3.  Dynamical  Systems  Explanations  of  Turbulence 

A  third,  and  somewhat  overlapping,  area  of  current  research  in  instability  and 
transition  is  the  application  of  dynamical  systems  theory  to  examine  the  temporal  behavior 
of  flows  in  transition.  These  theories  look  at  the  temporal  evolution  of  the  most  unstable 
eigenvalue  of  the  Navier-Stokes  equations,  linearized  about  the  current  state.  The 
bifurcation  of  the  flow  solution  from  one  state  to  another  can  be  explained  as  eigenvalues 
of  the  solution  crossing  from  the  stable  regime  to  the  unstable  regime. 

Currently,  several  scenarios  of  bifurcation  are  thought  to  pertain  to  the  transition 
of  flows  from  a  laminar  state  to  a  turbulent  state.  The  Ruelle-Takens-Newhouse  scenario 
occurs  as  a  sequence  of  Hopf  bifurcations  as  the  Reynolds  number  is  increased  slowly  in 
time.  That  is,  a  stable  flow  becomes  unstable  when  the  first  eigenvalue  pair  crosses  the 
growth-rate=0  axis.  The  flow  will  then  bifurcate  to  a  periodic  flow  with  a  single 
frequency.  This  process  is  repeated  when  a  second  eigenvalue  pair  becomes  unstable  and 
adds  a  second  frequency  to  the  flow,  making  it  quasi-periodic.  Dynamical  systems  theory 
suggests  that  when  a  third  frequency  appears,  incommensurate  with  either  of  the  first  two, 
the  nonlinear  interactions  will  quickly  fill  in  the  frequency  spectrum  and  a  non-periodic 
state  will  occur.  This  state  is  called  chaos.  A  second  scenario,  the  Feigenbaum  scenario 
explains  an  alternate  route  to  chaos  through  repeated  period-doubling  bifurcations.  This 
scenario  and  others  are  described  in  Guckenheimer  (1986). 

These  theories  can  be  investigate  and/or  used  to  explain  the  early  stages  of 
transition  in  channel  flows.  The  most  notable  use  of  the  above  explanations  is  in  Taylor- 
Couette  flow,  or  the  flow  between  two  concentric  cylinders  driven  by  the  motion  of  one  or 
both  of  the  cylinders.  This  is  not  a  channel  flow  due  the  difference  in  wall  boundary 
conditions,  but  could  certainly  be  investigated  using  the  methodology  developed  in  this 
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study  by  simply  applying  non-zero  wall  velocities  and  by  removing  the  pressure  gradient 
which  drives  channel  flows.  See  the  simulation  of  Streett  and  Hussaini  (1991). 
Applications  to  channel  flows  induced  by  a  pressure  gradient  have  also  been  investigated, 
most  notably  curved  channel  flow  (Finlay,  Keller  and  Ferziger,  1988).  Guzman  and  Amon 
(1994)  show  that  the  2-D  flow  in  their  periodic  converging-diverging  channel  follows  the 
Ruelle-Takens-Newhouse  scenario  as  the  Reynolds  number  is  increased. 

There  are  many  directions  to  pursue  to  investigate  the  applicability  of  dynamical 
systems  approaches  to  explanations  of  transition.  Again,  the  first  step  is  to  reproduce 
some  of  the  results  found  in  the  literature  showing  bifurcations  to  quasi-periodic  flows  or 
demonstrating  period  doubling.  The  next  step  is  to  search  for  new  flows  in  which  these 
phenomena  occur.  The  effect  of  spanwise  disturbances  on  the  bifurcation  process  should 
also  be  investigated  as  it  is  probable  that  certain  2-D  flows  would  never  arise  in  a  practical 
flow  due  to  3-D  instability. 

5.4.  Summary 

In  summary,  the  study  of  the  early  stages  of  transition  has  many  areas  which  need 
to  be  pursued.  The  general  code  described  in  this  study  can  be  used  for  a  wide  variety  of 
problems  which  can  address  these  areas.  The  next  phase  of  this  study  will  complete  the 
development  of  the  3-D  code,  and  then  will  select  one  of  these  areas  to  investigate  in 
greater  depth. 
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